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Parti.  Project  Summary 


This  report  summarizes  the  work  performed  under  DOD  SBIR  contract  #  DAAL06-92-C- 
0006,  to  develop  of  a  computational  algorithm  for  modeling  the  acousdc  scattering  from 
a  turbulence  distribution,  which  is  modeled  as  a  collection  of  discrete  turbules  with  a 
discrete  turbule  size  spectrum.  The  methods  adopted  here  build  on  the  prior  structure 
function  analysis  of  Tatarskii,  and  incorporate  this  approach  into  a  mean  field  method  for 
characterizing  coherent  and  incoherent  multiple  scattering  from  a  turbulence  volume. 
The  mean  field  treatment  is  an  extension  to  the  acoustic  case  of  prior  work  in 
electromagnetic  scattering  from  turbulence  by  one  of  the  authors  (Dr.  Resendes).^^  The 
object  of  the  current  effort  is  to  produce  a  phase  I  code  module  which  is  very  general,  and 
which  can,  under  phase  11  and  beyond,  be  adapted  to,  and  incorporated  in,  the  Army’s 
broader  acoustic  modeling  program,  with  particular  application  to  the  issue  of  scattering 
into  shadow  zones.  This  phase  I  SBIR  project  was  completed  on  November  18,  1992. 

Part  II  of  this  report  addresses  the  structure  function  approach  and  its  relationship  to  the 
computation  of  an  appropriate  eddy  number  density  distribution  as  a  function  of  eddy  size 
parameter.  The  work  discussed  here  builds  upon  prior  work  by  Dr.  Harry  Auverman‘1 
and  Dr.  George  Goedecke.5'7  Section  1  of  this  part  of  the  report  includes  a  discussion  of 
some  of  this  prior  work,  and  provides  a  background  for  the  balance  of  the  development  of 
the  structure  function  approach.  The  second  section  addresses  the  theoretical  derivation 
of  expressions  for  the  general  structure  function  and  for  the  power  spectral  density.  In 
section  3,  we  develop  a  refractive  index  structure  function  for  a  simple  turbulence  model, 
then,  in  section  4,  we  show  how  the  minimum  error  method  can  be  used  to  calculate  a 
number  density  dependence  on  eddy  size,  for  a  specific  eddy  model.  Section  5 
summarizes  the  results  of  part  n. 

In  part  HI  of  this  report,  we  begin  the  development  of  the  mean  field  approach  to  the 
multiple  acoustic  scattering  problem.  The  self  consistent  field  approach  is  applied  to  the 
multiple  scattering  of  acoustic  waves  from  a  turbulence  distribution  consisting  of  multiple 
scattering  sites  of  various  sizes  confined  to  a  turbulence  volume.  This  approach  provides 
a  general  and  systematic  procedure  for  carrying  out  a  perturbation  series  for  the  system  of 
scatterers  as  a  whole. 

The  averaging  process  for  determining  the  mean  field  is  model  dependent.  In  the  second 
section  of  part  HI,  we  present  the  turbulence  model  employed  for  these  phase  I 
calculations.  Subsequent  efforts  which  build  from  the  current  work  may  wish  to  employ 
a  different  or  more  complex  turbulence  model,  which  will  require  re-evaluation  of  the 
mean  fields. 
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After  describing  our  turbulence  model,  we  discuss  the  properties  of  a  single  scatterer. 
This  is  required  as  an  initial  step  in  the  construction  of  the  configuration  dependent 
equations  for  the  scattered  field  of  a  system  of  N  scatterers.  In  particular,  these  equations 
are  expressed  in  terms  of  the  exact  single  scatterer  result.  We  also  present  an  exact  result 
for  a  particular  single  scatterer  model.  Any  other  model  for  the  single  site  could  be 
employed,  and  approximate,  or  even  empirical  solutions  could  be  used  when  exact 

solutions  are  unavailable. 

Following  the  discussion  of  the  single  scattering  site,  we  proceed  to  the  consideration  of 
the  coupled  scattering  properties  of  a  collection  of  N  scatterers.  This  section  is  very 
general,  and  culminates  in  the  presentation  of  (1)  an  equation  for  the  total  wave  in  terms 
of  its  contributions  from  the  incident  wave  and  the  scattered  wave  from  each  of  the  N 
sites,  and  (2)  a  system  of  self  consistent  coupled  equations  satisfied  by  the  set  of  scattered 
waves  from  the  N  sites.  These  equations,  taken  together,  fully  describe  the  multiple 

scattering  problem. 


Following  the  development  of  the  equations  for  the  scattering  from  N  sites,  we  briefly 
discuss  the  probabilistic  and  statistical  aspects  involved  in  the  characterization  of  the 
configuration  of  the  N  scatterers.  We  then  proceed  to  the  evaluation  of  the  coherent  and 
the  incoherent  scattered  fields,  as  expressed  in  terms  of  the  mean  field.  We  next  present 
the  results  of  the  preceding  sections  in  a  plane  wave  basis.  Evaluation  of  the  acoustic 
signal  which  would  be  received  at  a  remote  detector  is  accomplished  by  means  of  the  far 
field  approximation  to  the  Kirchoff  integral. 

At  the  end  of  the  theoretical  development,  a  comparison  is  made  between  the  current 
formalism  and  the  "first  Bom"  treatment  of  Tatarski.  It  is  shown  that  Tatarski’s 
calculation  of  the  scattered  power  follows  from  the  more  general  treatment  employed 
here  by  keeping  lowest  order  terms  only,  and  making  further  simplifying  assumptions  in 
the  averaging  process. 

In  the  final  sections  of  the  effort,  we  present  the  numerical  arguments  employed  in  the 
computational  algorithm,  as  well  as  provide  the  text  of  the  code  itself.  We  also  discuss 
the  ways  in  which  this  code  can  be  incorporated  into  the  shadow  zone  analysis. 
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Eddy  Number  Density  Distribution  and  the  Structure  Function 
Introduction: 

A  direct  calculational  approach  for  the  eddy  number  density  distribution  is  developed, 
using  the  formalisms  of  the  index  of  refraction  structure  function,  and  the  minimum  error 
method  for  calculating  number  density  distribution  for  a  particular  model  of  turbulence. 
Much  of  the  work  presented  here  follows  closely  the  form  and  content  of  prior  analyses 
by  Dr.  G.  Goedecke,  from  which  it  in  part  derives  its  motivation.  Significant  departures 
are  made  from  Dr.  Goedecke's  work  in  that  assumptions  that  eddies  must  be  confined  to  a 
lattice-like  array  of  cells  which  subdivide  the  volume  are  dropped  in  favor  of  the 
presumption  that  any  given  eddy  may  be  found  anywhere  in  the  entire  turbulence  volume 
with  equal  probability.  We  also  depart  from  Dr.  Goedecke  by  choosing  to  adopt  an 
approximate  means  to  obtaining  the  number  density  distribution  from  a  model  which  has 
been  fitted  parametrically  to  the  theoretical  power  spectral  density,  rather  than 
constraining  ourselves  to  accepting  only  a  model  from  which  the  ideal  power  spectral 
density  itself  can  be  derived.  These  departures  permit  us  somewhat  greater  latitude  in  the 
development  of  an  application  oriented  algorithm  for  modeling  acoustic  scattering  effects 
for  a  turbulent  atmosphere.  We  begin  the  discussion  by  addressing  the  latter  departure. 


One  might  well  expect  to  be  able  to  derive  the  structure  function  and  the  corresponding 
power  spectral  density  of  the  fluctuating  refractive  index  in  terms  of  the  constituent 
eddies  which  make  up  the  turbulence  volume  which  exhibits  these  fluctuations.  There  is, 
however,  a  subtlety  here,  related  to  the  theoretical  status  of  the  structure  function  and  of 
the  turbulent  eddies.  While  at  least  in  the  inertial  subrange,  the  structure  function  is 
derived  from  physical  laws  and  principles,  a  superposition  of  eddies  of  differing  size  is  a 
highly  simplified  and  idealized  model  of  turbulence.  It  is  therefore  unlikely  that  the  true 
physical  structure  function  should  be  derivable  from  a  particular  model  collection  of 
eddies.  A  more  modest  expectation  might  be  that  an  appropriate  collection  of  eddies 
should  allow  one  to  approximate,  to  the  desired  accuracy,  the  physical  structure 
function. 


Mathematically,  this  observation  on  the  non-deriv ability  of  the  structure  function  from  a 
simplified  model  presents  itself  in  the  fact  that  it  is  not  possible  to  obtain  a  particular 
arbitrary  function,  such  as  f  ~  r^/^,  from  the  superposition  of  an  incomplete  set  of 
functions,  such  as  those  which  are  discussed  below  in  the  context  of  the  power  spectral 
density. 
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Lacking  closure  for  the  basis  functions,  the  Minimum  Error  Method  provides  a  means  for 
generating  an  approximate  representation  of  the  structure  function  in  terms  of  an  eddy 
superposition.  We  detail  below  a  method  for  obtaining  a  suitable  structure  function 
approximation.  The  general  outline  is  as  follows: 


We  begin  by  making  some  observations  on  prior  work.  Next  we  write  down  the  general 
theoretical  expressions  for  a  general  structure  function  and  a  general  power  spectral 
density.  This  is  followed  by  a  development  of  expressions  for  these  quantities  based  on 
an  eddy  model  of  turbulence.  The  Minimum  Error  procedure  is  then  used  to  determine 
the  optimal  eddy  number  density  distribution.  Finally,  we  present  a  brief  summary. 


Section  II.  1 

Observations  on  Prior  Work 


As  discussed  above,  the  work  conducted  in  this  program  is  rooted  in  the  prior  efforts  of 
Dr.  Harry  Auverman  and  Dr.  George  Goedecke.  This  prior  work,  both  at  and  sponsored 
by  the  Meteorology  and  Acoustics  branch  of  the  Army  Atmospheric  Sciences  Laboratory, 
emphasizes  the  development  of  a  reasonable  physical  characterization  of  acoustic 
scattering  from  turbulence,  using  a  model  built  from  the  superposition  of  constituent 
eddies.  During  the  course  of  review  of  this  earlier  work,  it  became  apparent  that  certain 
aspects  of  the  formalism  might  impose  more  constraints  than  were  completely  necessary 
for  the  development  of  an  acceptable  scattering  model.  This  section  addresses  these 
difficulties. 


The  general  formalism  which  was  earlier  developed  seeks  to  first  characterize  the  number 
density  of  constituent  eddies  as  a  function  of  eddy  size  and  then,  for  each  size,  constructs 
a  lattice  of  cells  within  the  turbulence  volume  such  that  each  cell  contains  exactly  one 
eddy  and  the  collection  of  cells  for  a  given  size  yields  the  correct  volumetric  eddy 
number  density.  This  procedure  ensures  that  the  correct  number  density  versus  eddy  size 
distribution  is  achieved  throughout  the  entire  turbulence  volume.  Each  eddy  is  presumed 
to  exhibit  an  index  of  refraction  fluctuation  which  is  localized  in  space  by  an  envelope 
function  which  decays  to  zero  in  some  manner  as  the  distance  from  the  eddy  center 
increases.  Examples  of  such  a  construct  would  include  a  Gaussian  envelope  function, 
which  decays  asymptotically  to  zero  at  large  radius,  or  a  "hard  sphere"  model  in  which 
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the  index  fluctuation  falls  abruptly  to  zero  at  some  particular  radius.  In  either  case,  the 
eddy  is  characterized  by  some  particular  size. 


In  addition  to  the  spatial  localization  of  the  refractive  index  fluctuation  about  a  centroid, 
the  cell  approach  involves  the  localization  of  the  centroid  of  the  eddy  within  the  cell 
volume.  While  this  aspect  of  the  model  leads  to  a  computationally  reasonable  description 
of  the  system  of  eddies  ,  it  appears  to  lack  physical  motivation,  since  it  constrains  the 
eddies  to  particular  regions  of  space  without  providing  any  particular  physical  force  to 
implement  this  constraint.  It  is  apparent,  for  example,  that  a  large  longitudinal  density 
fluctuation  in  the  overall  medium  could,  in  the  physical  world,  result  in  the  near  total 
depletion  of  eddies  of  a  given  size  over  a  scale  larger  than  the  cell  dimension.  Such 
depletion  is,  however,  specifically  denied  by  the  cell  formalism. 


In  addition  to  this  physical  difficulty,  there  is  also  a  computational  difficulty  with  the 
prior  treatment.  This  difficulty  arises  as  a  consequence  of  the  factorization  of  the 
averaging  procedure  over  the  turbule  distribution.  The  ability  to  factor  an  ensemble 
average  into  the  product  of  two  separate  averages  over  sub-ensembles  arises  in  the  event 
that  the  two  sub-ensembles  are  uncorrelated.  In  the  earlier  work,  an  initial  factorization  is 
done  into  separate  ensembles  for  eddies  of  differing  size.  A  subsequent  factorization  also 
separates  the  average  into  the  product  of  separate  averages  for  the  density  fluctuation  at 
different  sites  for  eddies  of  a  single  size. 

This  particular  factorization  is  a  relatively  commonplace  practice  in  the  computation  of 
multiple  scattering  effects,  and  is  done  in  an  effort  to  obtain  a  sufficiently  simple 
description  of  the  system  that  one  might  be  able  to  find  a  solution,  but  the  assumptions 
carry  with  them  a  secondary  consequence.  The  factorization  is  only  possible  if  the  eddies 
are  uncorrelated,  since  otherwise  the  variables  are  not  independent.  This  particular 
factorization  is  therefore  equivalent  to  stating  that  no  two  eddies  communicate  in  any 
way,  or  similarly,  that  the  eddies  are  entirely  uncorrelated.  If  the  eddies  are  uncorrelated, 
then  once  again  there  is  no  longer  any  particular  reason  to  believe  that  they  should  be 
confined  to  a  particular  region  of  space  surrounding  a  lattice  point  One  would  be  more 
inclined  to  believe  that  any  particular  eddy  has  an  equal  probability  of  being  found 
anywhere  within  the  entirety  of  the  turbulence  volume.  In  this  case,  the  probability 
function  which  describes  the  likelihood  of  the  eddy  being  at  a  particular  point  in  space 
would  be  constant  throughout  space.  The  corresponding  Fourier  transform  of  this  spatial 
probability  distribution  would  then  be  a  delta  function  in  k  space. 


One  final  observation  which  was  made  during  the  review  of  the  earlier  work  has  to  do 
with  the  small  r  dependence  of  the  derived  structure  function.  The  structure  function  can 
be  rather  generally  described  as  being  of  the  form 
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D/r)  =  I  dk  <I)(k)  (1  - 


A  derivation  of  this  form  and  a  discussion  of  the  assumptions  which  lead  to  this  particular 
form  are  presented  below.  This  function  has  a  quadratic  small  r  dependence,  provided 
only  that  there  is  a  maximum  value  of  k  for  which  0(k)  is  non  zero,  owing  to  the  small 
angle  properties  of  the  sine  function.  If  <I)(k)  is  zero  for  all  k  larger  than  about  l%li  , 
where  I  is  the  small  end  of  the  inertial  subrange,  then  all  sine  contributions  to  D  will  be  in 
the  small  argument  limit  for  r<f ,  and  the  function  will  exhibit  quadratic  behavior.  A 
consequence  of  this  observation  is  that  it  will  not  be  necessary  to  compute  a  particular 
<I>(k)  associated  with  the  small  r  region  of  D. 

Section  IL2 

Expressions  Derived  from  Theory  for  the  Structure  Function 
and  the  Power  Spectral  Density 


In  analyzing  the  spatial  structure  of  meteorological  fields,  it  is  appropriate  to  apply  the 
method  of  structure  functions.  The  concept  of  structure  function  arises  from  a 
consideration  of  a  random  function  whose  mean  value  changes  slowly  and  smoothly, 
rather  than  remaining  constant. 


A  difficulty  which  arises  in  the  analysis  of  such  functions  is  that  it  is  not  immediately 
obvious  where  to  define  the  boundary  between  true  variation  of  the  mean  value  and  very 
slow  fluctuations  about  that  mean.  Whenever  the  mean  value  <  f  >  changes,  we  can 

consider  instead  the  difference  F  =  f(r  )-f(r  +  ^).  For  values  of  ^  which  are  not 
too  large,  slow  changes  in  f  do  not  affect  this  difference  ( at  least  approximately  ). 

The  structure  function  is  defined  by 

DKti,t2)  S  (  [f(rl)-f(i1))2>  .  n.2.1 

The  difference  between  the  values  of  f(  r  )  at  two  points  Tj  and  r2  is  chiefly  affected 
only  by  inhomogeneities  of  the  field  with  dimensions  which  do  not  exceed  the 
distance  I  rj  -  I  .  If  this  distance  is  not  too  large,  then  the  largest  anisotropic 

inhomogeneities  have  no  effect  on  f(r^)-f(r2),  and  DfC?!,"?! )  can  depend 
only  on  (  ).  By  contrast,  the  correlation  function  ) .  defined  by 
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II.2.2 


s  <  [  m)  -  <  f(ri)>  ]  [  «n)  -  <  f(ii)>  ]  > . 

is  affected  by  inhomogeneities  of  all  scales.  In  addition  to  this  homogeneity  property  if 
the  random  field  is  also  isotropic,  then  the  structure  function  depends  only  on  I  Pj  -  r2  I. 


A  locally  homogeneous  random  field  f(  r  )  can  be  represented  by 

I  OO  I  CIO  I  oo 

f(t)  =  f(Tf)+  I  I  I  ■'^)dO(lt),  II.2.3 

.00  -00  »oo 


where 

dO(ki)  da)*(i?2)  =  6(  i?i  -  ki )  dki  dk": 


and  ^  0  is  the  spectral  density  of  f.  A  general  form  for  the  structure  function  can 
be  obtained  for  the  locally  homogeneous  random  field  by  substituting  the  general  form  of 
the  field  (equation  II.2.3)  into  the  definition  of  the  structure  function  (equation  II.2.1), 
and  evaluating  the  average,  making  use  of  the  relation  II.2.4.  The  result  is. 

Df(ti,t2)  =  Df{ti-f2,ii^)  =  Df(  r  ,1^ )  =  Df(r) 

=  <  [fit)  -  > 


=  <  [f(t)  -  f(Tf )]  [f(f )  -  fit)]*  > 


.00  -00  -00 


4-00  -H»  +00 

I  J 


.00  -00  -00 


^ )  d<t>(l5)) 


=  f  f  f dk*, <I.(kl) 

^co  -00  -00 


4co  +00  +00 

=  j  J  J  di?i  <l>(ki)  ( 1  -  cos  i?i  •  "r  )  n.2.5 

•00  -oo  *00 
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If,  in  addition,  f  is  locally  isotropic,  then  Df<1?)  =  Df(  r  )  and  <I>(k )  =  <I>(k),  and 
relation  II.2.5  can  be  further  simplified.  In  particular, 

+00  +00  +00 

Df{r)  -  2  j  J  J  dk  sin0  d0  d(j)  0(k)  [1  -  cos(kr  cos9)],  II.2.6 

•OO  -OO  >00 

where  0  is  the  angle  between  and  .  On  performing  the  integrations,  this  gives 


DKr)  =  8)c  J  ( 1  -  ^  )  <I>(k)  dk.  n.2.7 

0 

We  therefore  come  to  the  overall  conclusion  that  for  a  locally  homogeneous,  locally 
isotropic  random  field,  the  structure  function  is  of  the  form  II.2.7.  As  discussed  above, 
this  function  has  a  quadratic  small  r  dependence  if  <I>(k)  posesses  a  large  k  cut  off, 
independent  of  the  detailed  character  of  the  spectral  density  function  <l>  .  In  the  next 
section,  we  apply  the  formalism  of  this  section  to  an  eddy  model  of  turbulence,  to 
characterize  the  acoustic  response  of  the  system. 

Section  nj. 

Refractive  Index  Structure  Function  for  a  Simple  Turbulence  Model 

As  a  model  of  turbulence  we  consider  a  collection  of  localized  spherical  eddies  whose 
size  parameter  forms  a  discrete  spectrum.  It  is  assumed  that  this  description  can  be  made 
arbitrarily  close  to  a  continuum  size  distribution  by  increasing  the  number  of  sizes 
allowed  in  the  discrete  size  spectrum.  We  denote  the  number  of  sizes  as  Ng,  and  the 
number  of  eddies  of  each  size  by  N^.  The  distribution  of  eddies  is  assumed  to  be 
completely  random  within  a  cubical  volume  V  =  L^.  It  assumed  that  the  individual 
eddies  may  each  be  found  anywhere  in  the  volume  V. 

The  refractive  index  for  this  turbulence  model  can  be  taken  to  be  a  uniform  background 
value,  assumed  to  be  unity,  with  a  superposed  fluctuation  about  this  background  value 
associated  with  each  of  the  eddies.  We  assume  that  all  eddies  of  the  same  size  have  the 
same  functional  form  for  their  associated  localized  refractive  index  fluctuation.  For  the 
chosen  decomposition  of  the  eddy  population,  this  statement  can  be  written  as 

-y  Ns  Na  ^ 

n(  r )  -  1  =  X  X  fa(r-rn„),  n.3.1 

(X  =  1  =  1 
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where  n(  Ir  )  is  the  refractive  index  of  the  turbulent  system,  is  the  position  of  eddy 
number  na  of  size  parameter  a,  and  f^  is  the  function  which  describes  the  localized 
refractive  index  fluctuation  associated  with  an  eddy  of  size  parameter  a , 

The  correlation  function  for  this  refractive  index  model  is  given  by 

U(r  ,  r’)  S  <  (n(r)  -  1)  (n(  r')  -  l))  (generally) 


Ns  Ns  No  Nfl  ,  , 

=  111  1  (  fa(r-rno)f(3(r'-rmJ>. 

a=lP  =  lno=lmp  =  l  P 

The  structure  function  may  be  expressed  in  terms  of  the  correlation  function  as 

D„(r,r-)^  <[  (n( r ) - 1)  - (n(  r') - 1) f > 

=  11(7,-?)  +  -  2n(?,t') 


n.3.2 


11.3.3 


To  evaluate  the  structure  function  for  the  selected  eddy  model,  we  first  evaluate  the 
correlation  function.  It  is  convenient  to  separate  relation  II.3.2  into  two  terms.  The  first 
contribution  is  that  portion  which  is  diagonal  in  the  size  parameter,  i.e..  for  which  P  =  ct . 
The  second  contribution  is  from  the  off  diagonal  terms.  This  yields 


^(r  ,  t') 

Ns 

+  1 

a-i 


N, 


N, 


N, 


=  I  I  Z  <  fa(7  -  )  fa(7'- r-jn  J> 

0=1  na=l  010=1  “ 

Z  ^  <  fo(r  ■  iVio )  C)>  •  n, 

6  =  1  na  =  1  mp  =  1  P 

(Ma) 


3.4 


A  similar  separation  may  be  conducted  for  the  eddy  indices  n  and  m,  for  a  given  size 
parameter.  This  yields 
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II.3.5 


+ 


Ns 

I 

a  =  1 


N. 


(P5^a) 


At  this  point,  it  is  appropriate  to  begin  making  simplifying  assumptions.  We  will  follow 
here  the  same  approach  as  that  taken  by  Dr.  Goedecke  and  others.  The  first  assumption 
we  will  make  is  that  averages  over  turbules  of  differing  size  are  uncoirelated.  This  will 
allow  factorization  of  the  last  term  in  relation  II.3.5.  The  result  is 


Ns  N„ 


a*i  nrt  =  i  “  ® 


■a 


Ns  N„ 


Nr 


+  I  I  I 

a  =  1  n„  =  1  nv,  =  1 

(n^j  ^  mQ() 


<  fa(r -*na)fa(r'-rm^)> 


Ns  Ns  Na  1^ 

+  S  I  I 

a  =  1  b  =  1  Urt  =  1  mn  =  1 

(Na)  ^ 


<  fa(r -i^„)  ><f|3(r'-i^p)> 


— >1  — > 


We  next  make  the  further  simplifying  assumption  that  averages  over  different  sites,  even 
of  the  same  size,  are  uncorrelated.  This  gives 


Ns  Na 


^7 , 7’)  =2  2  <  fa(7  -  r„„ )  )) 


a  =  1  not  =  1 


Ns  Na 
+  11 


N, 


a 

I 


a  =  1  na  =  1  nia  =  1 

(met  net) 


<  fa(r -rJia)  >  <  ^aCr'- r^^)) 
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Ns  Ns 

a=i  B=l 

(Na) 


Na 

I 

Hot  -  1 


mp 


<  fa(r -rn^)  ><fp(r'-rmp)>. 


It  should  be  pointed  out  here  that  these  particular  assumptions  amount  to  saying  that  the 
likelihood  of  finding  a  particular  refractive  index  fluctuation  at  any  given  place  is 
completely  unrelated  to  the  likelihood  of  finding  any  other  such  fluctuation  at  any  other 
place.  This  situation  is  really  not  very  likely,  and  these  assumptions  are  in  all  probability 
violated.  Nevertheless,  we  make  these  assumptions  anyway,  in  an  effort  to  obtain  a  base 
solution  from  which  to  build  a  more  sophisticated  description. 


We  continue  with  the  development  of  the  equations,  by  adding  and  subtracting  like  terms 
so  as  to  complete  the  sums: 


\i(r  ,  r') 


Ns 

I 

a  =  l 


Na 

I 

"0  =  1 


<  fa(r 


No  Na  Na 

+  I  I  I  <  fa(r  "  ina )  >  <  ^aC r'-  ^J> 

a  =  lna=lma=l  “ 

-  I  I  <  fa(r  - r^a )>  < ^cxC r'- *n_)> 

a  =  1  na  =  1  “ 

Ns  Ns  Na  ^ 

+  1111  <  fa(r-4)  XfpCt’-fL)) 

a=l  p  =  l  na=l  mp  =  l  o  k  -p 
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Which  can  be  rewritten  as 


Ns 

-  1 

a=  1 


l'  [  S  <fa(l'-rlj)][  I  ,  <  fo(r'- rm„))] 

a  *  1  n^  =  1  ^  “ 


-  S,  I,  [<  fa(r- 4  »][(«■?'- 4))] 


a  =  l  na=l 


Next  we  simplify  the  notation  by  introducing  the  functions 
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which  yields 


H(r  ,?')=!  I  (  fa(r  -  K  )  fo(r'-  rj  )) 

a  =  i  11^  =  1  “  “ 

+  |.K(r)][s<,(?')] 

OC  ^  1 

-  I  I  <  fo(r  -  iv,  J>  < «  r'- ?„„)> 

a  =  1  tio  =  I  “ 

+  S(  r  )  S(r') 


which  simplifies  to 


+S(r)S(r')  . 


The  corresponding  results  for  p(  r,  r)  and  )i(r',r')  are  obtained  by  direct 
substitution,  and  the  results  entered  into  equation  11.3,3,  to  obtain  the  structure 

function  D(  Ir  ,1^') . 


In  order  to  evaluate  the  averages  in  the  above  equations,  it  is  convenient  to  introduce  a 
Fourier  series  representation.  The  plane  wave  basis  described  by  the  set  of  functions 


I 

is  a  complete  orthonormal  set  on  the  interval  [  j  ].  We  may  therefore  expand  any 
function  of  defined  in  the  volume  V  =  according  to 


,n=0,±l,  ±2,  ±3... 


n.3.7 


f(r) 


4.00  +00  4.00 

2  s  s 

k^=-oo  ky=-oo  kj=-oo 

y  ^  ■?) 

.1  fj  e 


i(t  ■?) 


n.3.8 


where 


~  1  f  ^  -i  r  •  r 

fj  =  v  J  fd?)  e  n.3.9 

V 


and 


-00  <  ^j<  00,  {integers} 


For  this  basis  set,  the  completeness  and  orthonormality  are  respectively: 
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II.3.10 


=  6(r  -  r') 


and 


-i(t.  t’)  "? 


=  ?  f.f  • 


II.3.11 


Making  use  of  these  last  few  results,  it  is  now  possible  to  evaluate  the  averages  in  the 
functions  0(1^,"?')  and  )J.(  r ,  r').  The  average  of  the  function  r  -  Tn^  )  is  given  by 

<  fa(  r  ■  f'na )  >  =  J  « r  -  %  >  Pa(^na )  • 

V 

where  )  is  the  probability  that  the  eddy  iIq^  is  located  at  the  position  fnct* 

Making  use  of  the  result  n.3.8,  this  can  be  rewritten  as 


<fa(r -*na)> 


=  1 


I 

t 


Ok 


i  1^- 


-  r 


Oa 


Pa(rJia> 


=x 

It 


i  k-  r 


ok 


J  dV: 


-4 

na 


•i  k 


'na 


Pa(*na)'  n.3.12 


At  this  point,  we  make  our  departure  from  Dr.  Goedecke's  treatment  We  will  make  the 
assumption  that  the  eddies  are  not  confined  to  localized  cells,  but  rather  are  free  to 
wander  anywhere  in  the  volume  V  with  equal  probability.  This  assumption  can  be 
expressed  as 


Substituting  this  into  the  integral,  we  obtain 


n.3.13 


<fa(r -rJ,o^)>  =E 


1  k-  r 
e 
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II.3.14 


=  |: 


=  f 


a-a 


Similarly,  we  find 


<  fo(7-4)  ><fa(r’-rm„)> 


=  <  X  f-.  e 


n.)  y  f 


K  t 


^6k 


•it*  X 

1  k  •  (  r  -  fn^  ) 


=  si;  f- 


t  t 


=  X  ft  f -t  e 

^  (xk  a  -k 


.  -4  .  -4  \ 

1  k- (  r  -  r  ) 


=  S  f  V.  f. 


•  r-^ 


II.3.15 


Substituting  11.3. 14  and  n.3.15  into  n.3.6,  we  find 


2?  Na  r  „  ^ 

|4(r.r')--  S  X  I,  X  ^  f -t 

a  =  inct  =  i  V.  “■* 


1  k  R  T  T  1 

®  ■  oS  oa-l 


Ns  Na 


+  S  S  f_T{ 


Ns  N| 


a=l  n„  =  l  P  =  1  np  =  l 


20 


II.3.16 


Nc 


=  X  N, 


a  =  1 


a 


I 


u  f 


i  k 


ok 


Similarly,  we  have 

Nj 

"r ,?)  =  ’r'.'r')  =  S  Na  Z 

a  =  1  k 


ai 


Ns 

X  Na 
a=i 


n.3.17 


Substitution  of  the  correlation  function  into  equation  II.3.3  for  the  structure  function, 
gives 


Dn(r.r')  = 


Ns 

-  2  X  Ntt 

a  =  i 


I 


ok  a 


1c 


U.3.18 


which  can  also  be  written  as 
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II.3.19 


N, 


Dn(  ■?,■?')  =  2  I  Nq  Z  ^  ^  ^  ^  ’ 


a=  1 


because  the  sum  is  over  all  kx,  ky,  k^  from  -o®  to  +««,  and  (  is  even. 


If  we  now  take  the  limit  as  L  ->  «>,  the  result  is 


(2Jt)3 


Jd’S. 


which  yields 


D„(  7.?')  =  li ,  Na  ^  1  ft  U  (  1  -  “S  S  )  '  11.3. 

Which  is  of  the  same  form  as  Tatarski  equation.  1.41: 

Dn(  r  .t')  =  2  J  d^ic  0(ii )  (  1  -  COS  ii  •  S  ), 
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with 


II.3.21 


section  II.4 


Minimum  Error  Method 


In  this  section  we  apply  the  minimum  error  method^’^  to  determine  the  optimal  number 
density  distribution  corresponding  to  a  specified  theoretical  power  spectral  density 
function.  The  total  squared  error  of  the  model  power  spectral  density  (k)n,o<iel  relative 
to  a  theoretical  power  spectral  density  <I>(k)theory  given  by 

kmax 

E(d>(k)n,odel )  “  J  [‘^0^)lheory  "  ^0^)modelP^  i 

kmin 


22 


In  relation  II.4.1,  <I>(k),},eory  is  derived  from  the  structure  function,  while  <I>(k)model 
be  expressed  as  a  sum  of  the  form 


Ns 

<I>(k)model  =  I .  %  .  n.4.2 

a  =  1 

where  the  n^  are  the  number  densities  of  eddies  of  Size  parameter  a,  which  serve  as 
expansion  coefficients  to  be  optimized.  The  are  particular  to  the  model,  and  will  be 
discussed  further  below. 


The  quantity  E(<l>(k)model)  is  minimized  if 


—  =0,  for  P  =  1,2,...,  Ns. 
dnp 


Performing  the  differentiations  in  II.4.3  yields 

8E 


0  —  ^  =  J  dk  2[<l>(k){},eQ^  “  ^(k)model]  ‘^0t)model 


P  k 


9np 


mm 


n.4.3 


On  evaluating  —  ^(k)niodel  »  die  result  is 


9n 


P 


?  Ns 

Sn„>P„(k)  =  4'p(k). 

dnp  dnp  a  =  1 


Substitution  of  this  result  for  the  derivative  yields 


'^max 


0  -  J  dk  [^(k)theory  “  ^(k)inodel]  , 
J^min 


which  on  substitution  from  n.4.2  yields 


Ng 


max 


I  no  j  dk  foOc)  >Fp(k)  =  J  dk<I>(k)n,«,,4'p(k) 

Ot  -  1  . 

Kn 


^min 


ki 

kinin 
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At  this  point  it  is  convenient  to  introduce  the  usual  definition  of  the  inner  product  as 


('Pa.'i'p)  =  1  dk  -f  „0c)  'Pp(k) , 


n.4.4 


so  that  the  previous  result  may  be  compactly  written  as 
Ns 

I  na  (  '{'a  .'J'p  )  =  (  .^^^theory)  .  for  p  =  1,  2, Ns-  n.4.5 
a  =  1 


The  Uniform  Sphere 


As  a  definite  example,  we  consider  the  case  of  the  homogeneous  uniform  sphere.  The 
refractive  index  fluctuation  is  given  by 


Aa.  >  r -*na 

0,  otherwise. 


n.4.6 


The  Fourier  transform  of  the  index  fluctuation  is  given  by 


e ^ ^  faC? ■  rl J 


=  ^  Jdiji/sinede  Jr2dr 


0  0 


=  ^  ^  [  sin  (kaa)  -  kaa  cos  (kaa)] 


n.4.7 


with  the  properties  that 


f  ^  =  f  -♦  =  f  ,  ,  f ,  is  real,  and  f .  =  f ,  . . 
ok  a|kl  ok’  ok  ’  ok  a(-k) 


Making  use  of  equation  n.3.21,  we  have 
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<D(k)  = 


V 


N, 


(2JI)3  a?  1  ■ 

=  ^  1  [  k3  ^  '*  Oc^a))]  ^ 


=  I  no  (sin  (kact)  -  kaa  cos  (kaa))^, 

a  =  l 


n.4.8 


where  use  has  been  made  of  the  relation 


H 


Y  — >  constant  s  as  V  . 


We  may  now  make  the  association  with  equation  II.4.2,  identifying  that 


'Fa  =  (sin  (kaa)  -  kaa  cos  (kaa))^. 


n.4.9 


The  relations  11.4.5  form  a  matrix  equation  with  the  ( 'Fa,'Fp)  as  matrix  coefficients 
relating  the  unknown  vector  to  the  vector  (  'Pp  »^theory)  •  choice  of  ^theory  be 
dictated  by  the  structure  function.  We  have,  for  example,  <I)  oe  k'“^  for  the  Kolmogorov 
spectrum  in  the  inertial  subrange. 

Within  the  confines  of  the  uniform  sphere  model,  the  choice  of  magnitude  of  refractive 
index  fluctuation,  Aa,  for  each  size  bin  is  an  unspecified  parameter.  It  is  still  possible, 
however,  to  solve  the  system  of  equations  for  the  products  of  the  number  densities,  n^, 
with  the  fluctuation  magnitudes,  Aa-  The  ^a  can  be  written  as  the  product  of  a  scalar 
constant  with  a  function  of  ka<x.  This  gives 

^a=  CagCkao), 


where 
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2  Ag^ 
K 


and 


With  these  substitutions,  equations  n.4.5  can  be  rewritten  as 
Ns 

^  ^  ^ct  (  8a  ’Sp  )  “  ^P  (  Sp  ’^theory)  »  P  ~  !»  2, ...,  Nj. 


or 


Ns 

^  [n(jCa  ](  g(x  »§p  )  ~  (  8p  >*^theory)  »  P  =  1»  2, ....  Nj.  11.4.10 


The  system  of  equations  11.4.10  can  be  solved  for  the  products  n^Ca ,  without  prior 
knowledge  of  the  manner  in  which  the  magnitude  of  the  index  fluctuation  varies  with 
size.  We  include  below  a  simple  subroutine  for  evaluating  the  required  inner  products  for 
the  homogeneous  uniform  sphere  model.  The  routine  is  written  in  standard  FORTRAN, 
and  so  should  readily  compile  on  virtually  any  computing  machine,  using  virtually  any 
compiler. 


Numerical  Algorithm  for  Number  Densities 
program  fintk 

dimension  a(20),fint(20,20),fint2(20) 
c  obtain  an  upper  bound  on  the  number  of  steps  to 
c  allow  in  the  integration  process. 

write(6,*)'  input  integral  step  limit' 
read(5,’'‘)n 
c 

c  open  up  the  input  and  output  files 
open(2,file='fintk2.inp') 
open(  1  ,file='fintk2.out') 
c 

c  set  the  chosen  number  of  eddy  sizes 
m  =  5 
c 
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c  initialize  required  variables: 
c 

c  irrational  number  PI 

pi  =  3.141592653 
c 

c  lower  limit  of  k  space  integration  is  2PI/L,  where 
c  L  is  the  upper  end  of  inertial  subrange  eddy 

c  size.  L  is  about  100  meters,  =10000  centimeters. 

fk  =  2.*pi/10000. 
c 

c  upper  limit  of  k  space  integration  is  2PI/lo,  where 
c  lo  is  the  lower  end  of  inertial  subrange  eddy 

c  size,  lo  is  about  1  centimeter 

uplimk  =  2.*pi  . 
c 

c  beginning  step  size  dk  is  one  tenth  of  smallest  k 
deltak  =  .2*pi/10000. 
c 

c  initialize  each  integral 
do  700  k=l,m 
do  600  j=l,m 
600  rint(k  j)  =  0. 

700  fint2(k)  =  0. 

c  note,  The  fint(k  J)  are  the  inner  products 

c  between  the  g  functions, 

c  The  fint2(k)  are  the  inner  products 

c  between  the  functions  g  and  phi. 

c 

c  Read  in  the  chosen  eddy  sizes 
read(2,*)(A(k)4c=l,m) 
c 
c 

C - 

c 

c 

c  Perform  the  integration 
dol000i=l,n 
do900k  =  l,m 
do  800j  =  k,m 

c  note,  the  matrix  is  symmetric,  so  do  upper  half  only 

800  fint(k  j)  =  fint(k  j)+deltak*g(fk*a(k))*g(fk*a(j)) 

900  fint2(k)  =  fint2(k)+deltak*g(fk*a(k))*phi(fk) 

c 

c  increment  k  by  the  differential  dk 
fk=fk+deltak 
c 

c  modify  step  size  to  get  faster  run  time 

if(deltak  .le.  0.01*fk)deltak=1.02*deltak 
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c 

c  exit  loop  on  reaching  upper  limit 
if(fk.gt.uplirnk)go  to  1050 
c 

c  put  a  few  output  statements  into  the  loop  to  allow 
c  monitoring  of  accumulating  integrals. 

if(i.eq.  10)write(l  ,20)fk,((fint(k,j)  j=l  ,fti)4c=l  ,m) 
(fint2(ii),ii=l,m) 

20  format(  1  x,’k=',f  1 5 .6,/,6(5e  1 2.5 ,/),//) 

if(i.eq.  100)write(  1 ,20)fk,((fint(k,j)J=l  ,m),k=l  ,m) 
(fint2(ii),ii=l,m) 

if(i.eq.3()0)write(  1 ,20)fk,((fint(k,j),j=  1  ,m),k=  1  ,m) 
(fint2(ii),ii=l,m) 

if(i.eq.5()0)write(  1 ,20)fk,((fint(k,j),j=l  ,m),k=l  ,m) 
(fint2(ii),ii=l,m) 

if(i.eq.8()0)write(l  ,20)fk,((fint(k,j)J=l  ,m),k=l  ,m) 
(fmt2(ii),ii=l,nri) 

if(i.eq.900)write(l,20)fk,((fint(kj),j=l,m),k=l,m) 

(fint2(ii),ii=l,ni) 

if(i.eq.  lC)00)wTite(  1 ,20)fk,((fint(k  j),j=l  ,m),k=l  ,m) 
(fmt2(ii),ii=l»m) 

if(i.eq.  10000)write(l  ,20)£k,((fint(k  j),j=l  ,m),k=l  ,m) 
(fint2(ii),ii=l,m) 

1000  continue 
c  end  of  integration  loop 


c  fill  in  the  lower  half  of  the  matrix 
1050  dol200k=l,m 

do  1100j=k+l,m 
1100  fintq4c)=fint(k,j) 

1200  continue 

c 

c  perform  output  of  final  result  to  output  file 
write(l,30)fk 

30  format(10x,'k=’^10.3) 

write(l,10)((fmt(k,j)J=l  ,m),k=l  ,m),(fint2(i),i=l,m) 
10  format(lx,5el2.5) 

c 

c  write  result  to  screen 

write(6,10)((fmt(kj),j=l.m),k=l,m),(fint2(i),i=l,m) 

c 

c  that’s  it... 

stop 

end 

function  g(x) 

g  =  (sin(x)  -  x*cos(x))**2/x**6 
end 
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function  phi(x) 
phi  =  x**(-ll./3.) 
end 


This  routine  provides  the  required  matrix  elements  for  evaluating  that  number  density 
distribution  as  a  function  of  size  parameter  within  the  uniform  sphere  model,  for  any 
particular  choice  of  size  bins,  which  gives  the  best  least  squared  error  to  the  Kolmogorov 
spectrum.  Any  other  choice  of  refractive  index  model  can  be  used  in  place  of  the  uniform 
sphere,  by  simply  replacing  the  functions  g  and  phi  with  their  corresponding  equivalents 


Section  II.5 

Summary  of  Computation  of  Number  Densities 

Our  procedure  for  obtaining  the  eddy  number  density  distribution  corresponding  to  the 
atmospheric  index  of  refraction  is  summarized  below. 

1.  Starting  from  the  structure  function  of  interest,  as  shown  in  figure  1, 


Figure  1.  Index  of  Refraction  Structure  Function 


construct  its  theoretical  spectral  density.  This  step  can  always  be  accomplished  because 
the  correlation  function  and  the  spectral  density  are  Fourier  transforms  of  each  other  and, 
by  equation  n.3.3,  the  structure  function  is  expressible  in  terms  of  correlation  functions. 
Since  the  structure  function  consists  of  three  distinct  regions,  one  method  is  to  obtain  a 
spectral  density  for  each  region,  and  impose  a  continuity  requirement 
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2.  Based  on  a  model  of  a  single  eddy,  such  as  the  uniform  sphere  or  an  exponentially 
decaying  spherical  distribution,  calculate  the  model  spectral  density. 

3.  Apply  the  minimum  error  method  to  obtain  the  optimal  number  density  size 
distribution,  i.e.  the  number  of  eddies  per  unit  volume  of  each  size,  which  best  matches 
the  model  to  the  theory  derived  from  the  structure  function. 
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ni.l.  Introduction 


Part  in 

Acoustic  Scattering  from  Turbulence 


The  self  consistent  field  approach  is  applied  to  the  multiple  scattering  of  acoustic  waves 
from  a  turbulence  distribution  consisting  of  multiple  scattering  sites  of  various  sizes 
confined  to  a  turbulence  volume.  The  self  consistent  field  method  has  previously  been 
successfully  applied  to  the  multiple  scattering  of  vector  electromagnetic  waves.  The 
present  work  extends  the  methodology  to  the  multiple  scattering  of  scalar  acoustic  waves. 

In  constructing  the  self  consistent  field  method,  it  is  assumed  that  the  wave  scattered  by 
each  scatterer  is  proportional  to  the  wave  incident  on  that  scatterer.  The  wave  incident  on 
that  single  scatterer  is  assumed  to  consist  of  the  combination  of  the  wave  incident  on  the 
system  of  scatterers  and  the  scattered  waves  emitted  from  all  other  scatterers.  Similarly, 
the  waves  emitted  by  all  of  the  other  scatters  are  influenced  in  part  by  the  emitted  wave 
from  the  particular  single  scatterer.  In  this  way,  both  the  wave  incident  on  the  single 
scatterer  and  the  wave  emitted  by  that  scatterer  include  the  effects  of  scattering  to  all 
orders.  This  approach  provides  a  general  and  systematic  procedure  for  carrying  out  a 
perturbation  series  for  the  system  of  scatterers  as  a  whole. 

A  key  factor  in  the  methodology  advanced  here  rests  in  the  computation  of  a  mean  field. 
This  involves  evaluating  averages  over  the  configuration  of  scatterers  before  computing 
the  scattered  wave,  rather  than  evaluating  the  scattered  wave  resulting  from  particular 
configurations  and  later  averaging  the  results.  The  averaging  process  for  determining  the 
mean  field  is  model  dependent.  In  the  succeeding  section  below,  we  present  the 
turbulence  model  employed  for  these  phase  I  calculations.  Subsequent  efforts  which 
build  from  the  current  work  may  wish  to  employ  a  different  or  more  complex  turbulence 
model,  which  will  require  reevaluation  of  the  mean  fields. 

After  describing  our  turbulence  model,  we  discuss  the  properties  of  a  single  scatterer. 
This  is  required  as  an  initial  step  in  the  construction  of  the  configuration  dependent 
equations  for  the  scattered  field  of  a  system  of  N  scatterers.  In  particular,  these  equations 
are  expressed  in  terms  of  the  exact  single  scatterer  result.  Here  too,  the  final  results  will 
depend  in  part  on  the  model  selected  for  the  scattering  properties  of  the  single  site.  In  the 
third  section  of  this  part  of  the  report,  we  present  an  exact  result  for  a  particular  model  of 
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the  single  scatterer.  Any  other  model  for  the  single  site  could  also  be  employed,  and 
approximate,  or  even  empirical  solutions  could  be  used  when  exact  solutions  are 
unavailable. 

Following  the  discussion  of  the  single  scattering  site;  we  proceed  to  the  consideration  of 
the  coupled  scattering  properties  of  a  collection  of  N  scatterers.  This  section  is  very 
general,  and  culminates  in  the  presentation  of  (1)  an  equation  for  the  total  wave  in  terms 
of  its  contributions  from  the  incident  wave  and  the  scattered  wave  from  each  of  the  N 
sites,  and  (2)  a  system  of  self  consistent  coupled  equations  satisfied  by  the  set  of  scattered 
waves  from  the  N  sites.  These  equations,  taken  together,  fully  describe  the  multiple 
scattering  problem. 

Following  the  development  of  the  equations  for  the  scattering  from  N  sites,  we  briefly 
discuss  the  probabilistic  and  statistical  aspects  involved  in  the  characterization  of  the 
configuration  of  the  N  scatterers.  In  particular,  we  discuss  the  extent  to  which  the 
probability  of  finding  a  particular  configuration  is  or  is  not  factorable  into  separate 
probabilities  for  finding  a  particular  scatterer  at  a  particular  given  position. 

We  next  proceed  to  the  evaluation  of  the  coherent  and  the  incoherent  scattered  field. 
Rather  than  solving  explicidy  for  the  scattered  fields  and  then  averaging,  the  averaging 
operation  is  taken  first  in  order  to  find  an  approximate  equation  obeyed  by  the  mean  field. 
The  mean  field  is  a  much  simpler  object,  depending  on  the  observation  point  only,  and 
not  on  the  positions  or  states  of  the  individual  scatterers.  The  coherent  and  incoherent 
scattering  are  then  expressed  in  terms  of  the  mean  field. 

In  preparation  for  the  implementation  of  the  approach  in  a  numerical  algorithm,  we  next 
present  the  results  of  the  preceding  sections  in  a  plane  wave  basis.  Evaluation  of  the 
acoustic  signal  which  would  be  received  at  a  remote  detector  is  accomplished  by  means 
of  the  far  field  approximation  to  the  Kirchoff  integral.  A  comparison  is  made  between 
the  current  formalism  and  the  "first  Bom"  treatment  of  Tatarski.  It  is  shown  that 
Tatarski's  calculation  of  the  scattered  power  follows  from  the  more  general  treatment 
employed  here  by  keeping  lowest  order  terms  only,  and  making  further  simplifying 
assumptions  in  the  averaging  process. 


32 


In  the  final  sections  of  the  effort,  we  present  the  numerical  arguments  employed  in  the 
computational  algorithm,  as  well  as  provide  the  text  of  the  code  itself.  We  also  discuss 
the  ways  in  which  this  code  can  be  incorporated  into  the  shadow  zone  analysis. 

ni.2.  Turbulence  Model 

Our  model  of  turbulence  consists  of  a  finite  concentration  of  homogeneous,  discrete, 
spherical  scatterers  of  refractive  index  n  and  radius  a.  The  scatterers  are  permitted  to  be 
of  various  sizes,  but  all  scatterers  of  a  given  size  are  identical.  To  simplify  the 
calculations,  all  scatterers  are  assumed  to  be  confined  within  a  rectangular  volume,  but 
distributed  in  a  random  way,  with  no  correlations  between  scatterer  locations.  The 
independent  parameters  of  the  model  are:  1)  the  number  density  of  scatterers  of  a  given 
size,  2)  the  refractive  index  of  each  scatterer,  and  3)  the  radius  of  each  scatterer.  These 
quantities  are  fixed  by  experimental  data  on  the  turbulent  atmosphere,  supplemented  by 
theoretical  information  such  as  that  provided  by  the  index  of  refraction  structure  function. 

nL3.  The  Single  Scatterer 

The  self  consistent  approach  to  the  multiple  scattering  problem  is  a  bootstrapping 
procedure  which  begins  with  a  solution  for  the  scattering  of  a  wave  from  a  single  site. 
For  the  current  effort,  we  will  work  from  an  exact  solution  for  single  site  scattering,  using 
the  homogeneous  sphere  model  discussed  above.  An  approximate  solution  using  some 
other  model  could  also  be  used,  provided  only  that  the  scattered  wave  from  a  particular 
site  is  linearly  related  to  the  incident  wave  on  that  same  site. 

The  linearity  assumption  is  easily  represented,  using  the  mathematical  formalism  of 
Hilbert  spaces.  The  linearity  is  preserved  whether  we  use  a  configuration  space 
representation  for  the  wave,  a  momentum  space  representation,  or  any  other 
representation  which  lends  itself  to  the  selected  model.  In  particular,  if  we  let  the 
incident  field  be  described  by  the  simbol  l^j),  then  the  scattered  field  can  be  described  as 

a  linear  function  of  the  incident  field: 

l7ts>  =  tl)Ci>.  ra.3.1 

In  its  most  general  form,  the  scattering  operator  t  depends  on  the  incident  and  scattered 
wave  vector,  the  material  properties  of  the  scatterer,  and  on  the  geometry  of  the  scatterer. 
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but  not  on  the  fields.  The  fields  ITCj)  and  l;ij)  are  Hilbert  space  vectors,  and  the  scattering 
response  iT  is  an  operator  in  this  space.  The  equation  III.3.1  is  an  abstract  and  formal 
representation  of  the  scattering  process.  In  the  next  paragraph  we  address  the  explicit 
coordinate  space  representation. 

The  basic  equation  for  sound  propagation  in  a  moving  medium  can  be  written  (see 
Tatarski,  equation  5.1)  in  the  form 

V2p  -  -^  (-  +  Uj— P  =  0 , 
c  3t  Dxj 

where  P  is  the  potential  of  the  sound  wave,  the  uj  are  the  components  of  the  velocity  of 
the  medium,  and  c  is  the  velocity  of  sound.  If  we  assume  a  harmonic  oscillation  for  the 
sound  wave,  then  P  can  be  written  as 

P  =  ite-“. 


On  inserting  this  form  for  P  into  equation  III.3.2,  and  keeping  only  terms  to  first  order  in 
the  dimensionless  fluid  velocity  u/c,  we  get  the  linear  result  (Tatarski,  5.5) 


V^rc  +  k^Ti  =  -2ik 


\^7t+  k2Y7t, 


in.3.3 


where  k  =  co/c  is  the  wave  vector  of  the  sound  wave.  We  note  that  equation  III.3.3  is 
linear.  We  further  note  that  the  properties  of  the  medium  enter  into  the  equation  via  the 
fluid  velocity  T?,  and  the  temperature  dependence  of  the  sound  velocity  c.  For  simplicity 
in  this  initial  development  of  the  model,  we  will  drop  the  velocity  dependent  term  in 
III.3.3,  and  consider  only  temperature  fluctuations.  Treating  the  remainder  of  the  right 
hand  side  of  the  equation  as  a  source  term,  equation  in.3.3  may  be  rewritten  as  an 
inhomogeneous  integral  equation: 


T^)  -  *0(1^) -if  J  JV- 


7t( 


ni.3.4 


where  7Cq(7^)  =  e  ‘  ^ 
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The  second  term  in  III.3.4  describes  the  scattered  wave  both  close  to  and  far  from  the 
scatterer.  For  a  homogeneous  spherical  scatterer  of  radius  a,  equation  IIL3.3  (or  4)  has  a 
well  known  exact  solution.  In  this  case  n  satisfies  a  homogeneous  equation, 

V2ii  +  k27t=0,  ni.3.5 


outside  the  scatterer,  and  a  similar  equation  inside  but  with  a  modified  wave  vector.  The 
solution  everywhere  is  obtained  by  demanding  continuity  at  the  spherical  boundary.  For 
scattering  by  a  sphere,  we  take  advantage  of  the  symmetry  by  choosing  functions  n  that 
satisfy  the  homogeneous  wave  equation  in  spherical  polar  coordinates  (r,  0,  <j)). 


The  eigenfunctions  7t(r,  0,  <ji)  =  R(r)  0(0)  <!>((}>)  that  are  single  valued  in  (|)  (except 
possibly  at  points  on  the  boundaries  between  regions  with  ±‘ssimilar  properties)  and  finite 
at  0  =  0  and  0  =  it  are  given  by 


_  m 

^  (even) 
fodd'y 


cos(m(})) 

sin(m(})) 


Pn"’(cos0)  Zn(kr), 


in.3.6 


where  Pn'"(cos0)  is  the  associated  Legendre  function  of  the  first  kind  of  degree  n  and 
order  m  (  n  =  0,  1,  2, ...,  m  =  0,  1,  2,...,  n),  and  Zn(kr)  is  any  of  the  four  spherical 

Bessel  functions: 

jn(P)  =  ^n+l/2^P^ 

yn(P)  =  ^n+l/2^P^ 

hn^^^(P)  =  j„(P)  +  iyn(P) 

=  j^Cp)  -  iyn(P)  • 

Here  J  and  Y  are  the  half  integral  Bessel  functions  of  the  first  and  second  kind, 
’  n+l/2  n+1/2 

jn  and  yn  are  the  spherical  Bessel  functions  ( with  uitroduced  for  convenience),  and 
hn^^^  and  hn^^^  are  the  spherical  Bessel  functions  of  the  third  kind  (also  known  as 
spherical  Hankel  functions). 
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Scattering  from  a  Uniform  Sphere 


In  this  section  we  present  a  rigorous  treatment  for  scattering  from  a  localized  refractive 
index  perturbation  which  can  be  modeled  as  a  uniform  constant  value,  different  from 
unity,  inside  a  sphere  of  radius  a,  with  a  refractive  Index  of  unity  everywhere  outside  the 
sphere.  The  geometry  of  this  configuration  is  illustrated  in  figure  2. 


Figure  2.  Uniform  Sphere  Model 

Region  I  is  taken  to  be  the  region  of  unperturbed  atmosphere  outsied  the  sphere,  where 
the  refractive  index  is  unity.  Region  11  is  presumed  to  also  be  characterized  by  a  constant 
refractive  index,  but  which  differs  from  unity.  In  either  case,  the  propagation  of  sound  is 
characterized  by  a  source  free  (homogeneous)  differential  equation: 

Region  I:  (V2  +  k2)  Jtj  =  0 

Region  II:  (V^  +  q2)  7C2  =  0  in.3.8 


We  require  that  the  solution  should  behave  as  an  outgoing  spherical  wave  at  infinity,  in 
region  I,  and  that  the  solution  should  be  finite  at  the  origin,  in  region  H, 

At  the  interface,  the  two  regions  are  connected  by  appropriate  continuity  conditions.  In 
particular,  if  the  pressure  on  the  inside  of  the  interface  differs  from  that  on  the  outside, 
then  the  arbitrarily  abrupt  interface  would  be  subjected  to  an  infinite  acceleration.  We 
therefore  require  that  at  r  =  a,  the  pressure  must  be  continuous,  i.e.,  at  r  =  a,  pj  =  pjj. 
The  pressure  continuity  requirement  is  in  essence  equivalent  to  conservation  of 
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momentum;  if  the  fluid  were  to  flow  with  respect  to  the  turbule,  we  would  require 
conservation  at  the  boundary  of  the  total  momentum  p  +  pu^  (in  a  frame  fixed  to  the 


sphere).  The  second  requirement  is  that  the  normal  component  of  velocity  should  be 
continuous.  This  requirement  is  essentially  that  of  conservation  of  mass,  with  the 

generalization  for  the  flowing  fluid  pi?  +  ^  =  0,  since  we  do  not  have  sources  or 


sinks  at  the  interface. 


An  acceptable  interior  solution  can  be  constructed  from  those  polar  coordinate  solutions 
to  the  homogeneous  differential  equation  which  are  finite  at  the  origin: 


oo  n  \  m  T 

Jt2(r,e,(|))  =^)n(qr)  [an^PnCcosG)  (cos0)]. 


m.3.9 


Acceptable  exterior  solutions  are  constructed  from  those  solutions  which  behave  as  an 
outgoing  spherical  wave  at  infinity: 


;ci(r,e,(l))  hn^^\kr)[cn^Pn(cos0)+^Ij(cnniCOS  m(J)-KinmSin  m(t))  Pn'"(cos0)] 


ra.3.10 


In  addition  to  the  scattered  wave  ,  region  I  also  includes  the  incident  wave: 


JtincM’*!')  =  =  n5)  (2n+l)  i"  jnOo-)  Pn(cos0). 


m.3.11 


The  condition  of  continuity  of  pressure  at  the  interface  can  therefore  be  written  as 
-icop’  7C2  =  -icop(7ti+Jtinc)  •••  ( evaluated  at  r  =  a) 

or  equivalently, 


37 


II  m 

P’  n5  [anoPn(cos0)  +^Zj(anm  COS  m(|)  +  bnm  sin  m(j))P„  (cos0)] 

=  p  hn  Xka)  [cn^Pn(cos0)  +^Zj  (cnm  cos  +  dnm  sin  m(j>)  Pn  (cos0)] 

+  (2n+l)i"jn(ka)Pn(cos0), 

from  which  we  find,  on  equating  the  coefficients, 

p’  jnCqa)  ap^  =  p  [hn^^Wcn^+  (2n+l)  i"  jn(ka)  ] 

P’  jn(qa)  anm  =  P  hn^^\ka)  Cnm  ...  (m>l) 
p'  jn(qa)  bnm  =  P  hn  ^\ka)  dnm  ...  (ni>l) . 

Similarly,  the  continuity  of  the  normal  component  of  ii  can  be  written  as 

.  (evaluated  at  r=a), 

ar  ar 

from  which  we  find 

q  jn(qa)  an^  =  k  [h;^^\ka)  Cn^+  (2n+l)  i"  j;(ka)  ] 

qjn(qa)  anm  =  k'hn  Cnm  ...  (m^l) 
q  jn(qa)  f^nm  “  O^a)  dnm  ...  (mSl) ...  in.3.14 

For  the  general  k,  q,  and  a,  these  equations  make  contrary  demands  for  the  relationship 
between  anm  and  Cnm .  which  can  only  be  simultaneously  satisfied  if  anm  =  Cnm  =  0.  The 
same  argument  applies  to  bnm  and  dnm  .  We  therefore  conclude  that  the  7t(r,0,(j))  must  be 
of  the  form 


ni.3.12 


ni.3.13 
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r 


with  an^ 


7r(r,e,<|))=  i 


„5)  jnCqOan^PnCcose)  ;  r<a 


L  n5)hn^'W)Cn^Pn(cos0)  ;r^a, 
and  given  by  the  solution  to  the  coupled  equations 

P’  jn(qa)  an^  •  phn^^W)Cno  =  P  (2n+l)  i"  jn(ka) 
qjn(qa)  an^  -  k  h^^^^ka)  =  k  (2n+l)  i"  jn(ka) . 


IIL3.15 


ni.3.16 


Hilbert  Space  Formulation 

To  conclude  this  section,  we  return  to  the  fundamental  acoustic  wave  eauation  with  which 
we  opened,  but  this  time  continue  the  development  in  our  abstract  Hilbert  space 
formalism.  Using  this  approach,  we  will  rewrite  the  equation  for  sound  propagation  in 
the  absence  of  fluid  motion,  and  solve  it  exactly  ,  in  a  formal  manner.  Reproducing 
equation  131.3.3  (but  neglecting  the  velocity  contribution),  we  have 

(V2  +  k2)jt  =  k2Yjt, 

The  propagator  in  the  the  unperturbed  atmosphere  is  given  by 

(V2  +  k2)go  =  1 

Where  1  indicates  the  unit  matrix  or  identity  operator.  On  solving  this  equation  for  the 
free  propagator  go,  we  get 


in.3.17 


in.3.18 


^  =  (V2+k2r'. 


ni.3.19 


Equation  in.3.17  may  now  be  rewritten  in  Hilbert  Space  notation  as 

lo'^  )  =  k2  Y  lit  > . 

The  general  solution  to  in.3.20  is  given  by 

IJI )  =  IJCo  >  +  to  k2  ^  IJt  > , 


m.3.20 


m.3.21 
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where  IJto  )  satisfies  a  homogeneous  equation  which  can  be  written  in  the  coordinate 

space  representation  as  (V2  +  k2)7Co  =  0.  Equation  III.3.21  is  an  integral  equation 

t* 

for  Itc  )  ,  which  can  be  solved  by  iteration.  The  operator  -  is  responsible  for 

scattering,  and  is  therefore  a  scattering  potential;  for  simplicity,  we  rename  it  to  be  O’. 
Iterating  on  equation  in.3.21,  we  find 

>  +  go  v'  IJt  > 

=  )  +  go  v'  [  )  +  go  v'  [  )  +  -  ]] 

=  iTto  >  +  to  v’  iTto  )  +  (to  f  '^0  >  +  - 

=  Ij^o  )  +  to(  v'  +  O'  O'  + ...  )  Ijto  > 

=  teo>+toO'(l  ■  So«')''l’to> 


where 

=  IJto  >  +  go  S  > 

ni.3.22 

S^(l-goO’)'‘ 

ni.3.22 

=  0'+  O'tov’  +  v’to  O'toO'  +...  . 

in.3.23 

The  result  111.3.22  is  an  exact  solution  to  III.3.21,  valid  to  all  orders  in  the  scattering 
potential.  If  we  take  the  incident  field  to  be  given  by  ItIq  ) ,  then  the  term  go§  ItCq  )  is  the 
exact  solution  for  the  scattered  wave  and  is  therefore  equal  to  )  in  equation  III.3.1, 
provided  we  take  the  incident  field  to  be  lit, )  =  litp ) .  In  practice,  this  incident  wave  is 
generally  chosen  to  be  a  plane  wave.  To  summarize  these  results,  we  can  make  the 
identification 

t=  toS  =  tov'(l  - 

Successive  orders  of  approximation  to  are  given  by 

gov’  teo)»  etc. 

The  contribution  Itc  is  the  quantity  identified  as  the  Bom  approximation. 
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Section  IIL4 


The  System  of  N  Scatterers 

The  cornerstone  of  the  development  of  the  scattering  formalism  for  a  system  of  scatterers 
is  the  solution  to  the  single  site  scattering  problem  given  in  section  in.3: 

iTCs)  =  t  iTCj) .  ni.4.1 


When  considering  a  system  of  N  scatterers,  however,  the  wave  incident  on  a  given 
scattering  center  is  not  merely  the  wave  incident  on  the  system,  but  is  rather  the 
combination  of  that  incident  wave  and  all  the  scattered  waves  emanating  from  the  other 
N-1  scattering  sites.  Let  Ipj)  represent  the  acoustic  wave  incident  on  the  system  of 

scatterers,  and  let  l7ts(r^,'^j  )>  represent  the  wave  scattered  by  the  scattering  center 
located  at  the  site  ,  where  is  a  generalized  coordinate  comprising  both  the  position 
and  size  a.  The  total  acoustic  wave  can  be  represented  by  an  expression  of  the  form 

=  'Pi)  f  ))  ’  ^A.2 


which  explicitly  states  that  the  total  wave  is  the  sum  of  the  incident  wave  and  the  waves 
scattered  by  each  of  the  N  scattering  sites. 


The  total  wave  is  clearly  configuration  dependent  If  we  let  the  notation  }  denote  the 
set  of  scattering  center  locations  associated  with  a  particular  configuration,  then  this 
configuration  dependence  can  be  made  explicit  by  writing 

teTC^;(tj))>  =  ipi>+X  I  *,(■?,  C^j)'  )  >, 

where  the  notation  f^j)'  s  {  ,1^,^  allows  us  to  specifically 

extract  the  dependence  on  site  j,  and  emphasize  the  fact  that  the  wave  scattered  at  site  j 
depends  on  the  specific  configuration  of  the  other  N-1  scatterers. 
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The  wave  incident  on  scatterer  j  is  given  by 


=  !),,<■?)>- nw-: 

Making  use  of  the  single  scatterer  result,  IIL4. 1,  the  scattered  wave  from  site  j  can  be 
written  as 

))  =  t  )  litiC?)> .  ni  '*' 

If  we  substitute  the  result  in.4.4  into  the  relation  ni.4.2  for  the  total  wave,  we  get  the 
result 


=  IPi>  +  f 

If  we  further  substitute  III.4.4  into  in.4.3,  we  obtain 
system  of  linear  integral  equations 

l;c^(l^)>  =  Ipi)  +  iTt^Ci^)). 


|jt^(7^)) .  ni.4.5 

the  result  that  the  Itc^CT^))  obey  a 


Equations  in.4.5  and  in.4.6  constitute  a  set  of  self-consistent  coupled  equations  which 
completely  determine  the  multiple  scattering  problem. 


Section  in.5 
Statistics 

Let  pfC^  })  denote  the  configuration  probability  density  and  {d^bj}  the  volume  element 
for  the  N  scatterer  system.  The  probability  that  scatterer  1  will  be  found  in  d^bj 

about  ,  and  that  scatterer  2  will  be  found  in  d^b2  about  B2 »  ^nti  that  scatterer  3  will 


42 


be  found  in  d3b3  aboutB^ ,  and  ....  and  that  scatterer  N  will  be  found  in  d^b^  about 
is  given  by 

pCC^j })  {d3bj} , 

a  quantity  whose  integral  is  normalized  to  unity.  The  probability  density  for  a  particular 
subset  n  of  the  N  scatterers,  such  as  {Bi  ,^2  ,^3  ,..X :  n<N},  can  be  found  by  integrating 
the  probability  density  for  the  N  scatterer  system  over  the  remaining  variables  not  in  the 
particular  subset: 

p((  E5  ,B|  ,B|  ;  n<N))  =  J  d3b„..,  ...  d^bn  p((  )) .  ni.5.2 


We  will  denote  a  completely  random  distribution,  i.e.  one  which  factors  completely,  by 
the  probability  density 

Pn({1^))  =  Pi(^i)  P2(^)  PaC^)  -  Pn(^) 

If  the  probability  density  is  not  completely  random,  and  hence  does  not  factor  completely, 
but  rather  factors  only  partially,  this  can  be  represented  as 

where  the  last  factor  denotes  the  distribution  for  particles  n+1, ...,  N,  when  the  values  of 
{^ }  are  known. 


The  configuration  average  of  the  total  wave  is  given  by 


=  1  {d^bj}  p({1^  })  {"^  })>  • 

If  n  of  the  scatterers  are  held  fixed,  this  will  be  denoted  by  a  subscript: 


ni.5.5 
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=  J  Pn.„(  (S’n  )  mV;  (tj  ))>  . 


ni.5.6 


Ijit(  r )) 


{n} 


The  principal  use  which  we  will  make  of  expression  III.5.6  is  to  evaluate  the 
configuration  average  with  one  of  the  scatterers,  say  scatterer  k,  held  fixed,  in  which 
particular  case  expression  III.5.6  would  be  written  as 


d3bid3b2  ...d^bij.id^bic+i  ...d3bNPN.i(^/(~^k}')  })) 


ni.5.7 


Note  that  probability  distributions  may  be  converted  to  density  distributions  by 
multiplying  by  the  appropriate  power  of  the  number  of  scatterers,  for  example: 


n(bi )  =  N  pi(Bi );  nC^i  ,B2 )  =  p2(Bi  ,52 )  • 


Section  nL6 

The  Coherent  and  Incoherent  Waves 

It  is  generally  not  possible  to  obtain  the  total  wave  for  a  given  configuration  of  scatterers, 
nor  is  it  necessary,  since  we  are  interested  only  in  the  first  and  second  moments  of  the  N 
scatterer  probability  distribution.  In  this  section  we  will  use  the  statistical  relations  of  the 
preceding  section  to  perform  the  required  averages,  and  obtain  averaged  expressions  for 
the  coherent  and  incoherent  contributions  to  the  scattered  wave. 

CQhcr^nt  Scattering 

The  average  over  in.4.5  may  be  performed  using  relations  III.5.5  and  in.5.6,  along  with 
relation  ni.5.4,  with  n  =  1. 


Itct)  =  Ipi)  +  2  /  d^bj  piC'^j )  t  C^j ) 


where  the  object 
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l7C^)j  =  iTc’Cl^.'^j  ))j 

=  J  d3bid3b2...d3bj.id3bj+i...d3bN  PN-l('^j/ )’)  I  1)>  •  ni.6.2 

defines  the  effective  field,  and  represents  the  field  incident  on  the  particular  scatterer  j 
whose  position  is  known,  averaged  over  the  distribution  of  all  other  scatterers.  Note  that 

the  7^  dependence  is  still  implied,  but  has  been  dropped  from  the  explicit  notation  in 
relation  ni.6.1. 

The  effective  field  differs  from  the  total  field  by  the  field  emitted  by  one  scatterer.  If  the 
number  of  scatterers,  N,  is  large,  then  the  effective  field  is  approximately  equal  to  the 
coherent  field; 


=  iTtTCi^)) 


ni.6.3 


This  last  statement  is  true  for  each  of  the  N  scatterers,  so  that  when  this  approximation  is 
inserted  into  relation  III.6.1,  the  summation  over  j  is  simply  replaced  by  a  multiplicative 
factor  N.  In  particular  we  have 


Ijct)  =  Ipi)  +  t  Ijct>  » 


ni.6.4 


where 

t  =  N  J  d3bj  PiC^)  t(l^),  and  NpiC^j)  =  nC^). 


This  approximation  replaces  solving  the  system  of  N  coupled  linear  equations  for  the 
total  field  (relations  in.4.5  and  III.4.6)  by  a  single  integral  equation  for  the  coherent  field. 
Iterating  in.6.4  we  get 
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\TZj)  =  (1+  t  +  (  t  )2  +...)lPi> 

=(  1  -  t  ipi) .  ni.6.6 


Incoherent  Scattering 

The  density  and  directionality  of  incoherent  radiadon  is  given  by  the  <7^  L.l  r^ )  matrix 
elements  of  the  quantity  I  Kj  )( Jt-r  I.  We  now  proceed  to  obtain  an  equation  for  this 
quantity,  starting  from  the  basic  relation  IIL4.5  and  its  dual.  In  the  following,  we  will 
drop  the  explicit  7^  and  (Vj  }  dependences  from  the  explicit  notation,  leaving  those 
arguments  as  merely  implied.  Specifically,  we  get 

l7tT>  =  IPi)  +  f  tjlrc*)  in.6.7 

and  (kjI  =  <Pil  +  E  </IT/  .  ni.6.8 

On  taking  the  difference  between  the  average  of  the  product  and  the  product  of  the 
averages  of  equations  III.6.7  and  in.6.8,  we  get 


I  Jl-p  X  ^  ~  ^  ^  ^  "tj  I )  {  It  llj-  •  1 7C^ )  {it  l*fjj 

The  averages  of  the  diagonal  and  off  diagonal  terms  of  the  second  quantity  on  the  right 
hand  side  of  relation  in.6.9  are: 


tjiii'Xit'it'/  =/  d3bj  pi(Vj)  t(Vj)  ( Iitixit*!  )j  V(Vj) 


m.6.10 


=  Jd3bjd3bkP2(Vj,i^)tc^jX  “ 


4« 


Relations  III.6.9  through  III.6.11  are  exact  To  simplify  the  problem  we  make  similar 
approximations  to  those  made  in  the  coherent  case: 


(  Itc*)  )j  =  iTtx) 


(  I  TC*  >  <  JT*  I  )j  =  I  Kt  >  <  JCx  I , 

and  — - IT  -  ni.6.12 

I  )  ( «  I  )jk  -  I  JtT  )  ( I 

At  this  point,  it  is  convenient  to  introduce  three  super-operators,  R,  M,  and  Q  .  These 
super-operoperators  are  defined  for  and  operate  on  any  given  quantity  X  which  is 
independent  of  the  configuration  of  the  scatterers  and  their  states.  The  operators  are 
defined  by: 


R  X  =  t  X  t  ^ 

MX  =  N  J  d3bj  piC^j)  tC^j)xTWj) 
and  Q  X  s  N  (N- 1)  J  d^bj  d^bj,  PzCl^j  Mi  )t  C^j )  X  (t^(5j).  1^.6. 1 3 

Using  relations  in.6.12  and  the  super-operators  in.6.13,  it  is  possible  to  rewrite  equation 
in.6.9  as 

-  — -  -  — T  ■; —  in.6.14 

I  Jlx  X  ^  =  I  t^x  )  (  ^  Q  )  I  X  ^  -  R  I  JtX  )  (  ^ ' 

If  we  add  to  both  sides  the  quantity  R  ( I  Jtx  )  ( JCx '  *  1  JtT  X  *  )  >  we  get 

(  1  -  R  )  (  I  TCx  X  ^  ■  I  )  I  X  ^  » 
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which  is  readily  written  as  an  inhomogeneous  integral  equation; 


I JI7 )  ( I  =  I  JtT  )  ( ^  ^  I  Ttj )  ( JtT  I  , 


ni.6.16 


where 


L  =  ( 1  -R  )■*  (M  +Q  -R  ) .  ni.6.17 

Upon  iterative  expansion  of  ni.6.16,  we  obtain  the  following  power  series  for  the 
incoherent  radiation: 


iTtiXtil'  =(1+^-  +L'^  +  L'i  +...)  I  Ttj >  <  rtj  I 


=  (  1  -  L  )'^  I  JtT )  (  I  • 


From  III.6.18,  it  is  evident  that  the  coherent  radiation  term  Ijtj)  (Jtjl  acts  as  a  source 

term  for  the  incoherent  radiation  I  jcj  )  <  Jtj  I  •  The  composite  superoperator  L  generates 
successive  powers  of  incoherent  scattering.  On  neglecting  L  ,  or  equivalently  on 
considering  only  the  L  0  term  in  in.6. 1 8,  the  radiation  is  all  coherent  and  we  have 


.  ,  ,  .  V  /  .  Ul.0. 

I  Jlj  )  (  JCj  I  —  I  )  \  '  • 

The  operator!  generates  only  incoherent  scattering,  with  successive  powers  representing 
radiation  that  has  been  incoherently  scattered  once,  twice,  and  so  on.  The  M  term 
represents  the  purely  incoherent  addition  of  the  the  intensities  of  individual  scatterers, 
while  the  term  Q’R  represents  scattering  due  to  the  correlation  p2(15},b2)  * 
pl(Bi  )pi(B2  )  between  scatterer  locations. 
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Section  HI.? 


Coherent  and  Incoherent  Scattering  in  a  Plane  Wave  Representation 


In  this  section  we  recast  the  principal  results  in.6.6,  III.6.18,  and  in.6.19  in  a  plane  wave 
representation.  The  motivation  for  employing  plane  waves  for  the  system  of  N  scatterers 
is  that  the  configuration  dependence  appears,  in  this  representation,  as  a  simple 
multiplicative  phase  factor.  This  simplifies  the  averaging  procedure. 


The  Plane  Wave  Basis 

We  will  consider  a  cubical  box  V  =  ,  which  encloses  the  distribution  of  scatterers. 

Plane  waves  which  satisfy  periodic  boundary  conditions  at  the  walls  of  the  box  will  have 

— > 

wave  vectors  k  = — w — ,  with  o  =  cTj^  +  02^+035  •  ^2,  ^3  fill  integers.  We  can 

denote  a  complete  set  of  orthonormal  state  vectors  indexed  by  the  wave  vectors,  using  the 
Dirac  bracket  notation,  as  I  f  >  .  This  set  of  plane  wave  states  forms  a  complete 
orthonormal  set,  with  the  properties  of  completeness  and  orthnormality  given  by 

lik„)<i?,i  =  l 

a  O''  o 


and 


C>=  s, 


oo 


in.7.2 


If  we  introduce  a  momentum  operator  p  =  -i  ^  ,  then  this  operator  takes  on  wave 
vector  eigenvalues  so  that  (neglecting  the  size  parameter  ct  for  simplicity) 


in.7.3 


where  Tj  is  any  of  the  scatterer  locations.  We  use  this  property  to  compute  the  matrix 
elements  of  )  between  plane  wave  states.  In  particular,  we  have 

SO  that  the  matrix  elements  of  the  transition  operator  are 
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ni.7.5 


_  g  i  (ko  -  ko’)  •  rj  ^  I  t( "?  )  I  )  . 


The  meaning  of  IIL7.5  is  that  the  matrix  element  of  the  scattering  amplitude  for  a  single 

— ^  ^  ^ 

scatterer  between  plane  wave  states  I  k  ^ )  and  I  k  ^, )  for  a  scatterer  at  Tj  differs 

from  that  of  a  scatterer  at  the  origin  only  by  the  phase  factor  e  ^  ”  oV  *  *j  .  An 

important  point  here  is  that  the  random  variable  fj  appears  only  in  this  phase  factor. 


Coherent  and  Incoherent  Scattering 

We  will  begin  the  development  of  the  formalism  in  the  plane  wave  basis  by  computing 
the  expression  III. 6. 6  for  the  mean  total  coherent  wave.  This  is  accomplished  by 
computing  successive  terms  in  the  series  expansion  of  Itcj)  ,  and  performing  the 
infinite  sum.  We  will  assume  that  the  scatterers  are  distributed  completely  randomly 
throughout  the  volume  V  =  U  which  encloses  the  scatterers,  so  that  the  number  density 
n(rj)  is  given  by  n(rj)  =  Np(rj)  =  NA^ . 

Specifically,  equation  in.6.6  says: 


litj)  =  (l+  ^  +  (^)^+... )  Ipi) 


=(i-  t  )-‘ipi>. 


In  order  to  convert  this  to  a  plane  wave  basis,  we  can  invoke  the  definition  of  the  mean 
value  of  the  single  scatterer  transition  operator  to  obtain 


t  =N  xli^d3rj^(:?) 
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=  Ni;<f„i  “ho )  i?„)i  k„><f„i 


ni.7.6 


H  NZ  T^laXol  ni.7.7 

o  a  ' ' 

Making  use  of  this  result,  we  can  compute  the  second  power  of  the  average  transition 
operator: 

(T)2  =  T  •  T  =  (  NXT^IoXol  )x  (NjT^IrfXo'l  ) 

=  N2  2,  T,T^  loXoltfXo'l 
=  n2Z.T„T„,Io>5„^(o'I 
=  Z  (NT„)2  laXol 

The  same  argument  can  be  made  to  show  that  the  general  term  is 

Z(NT>la)(al, 

O  ^ 

so  that  the  full  series  (  1  -  1"  )’Ms  given  by 

(1-T)*l=  Z  [1+NT^  +  (NT^)2+...]  laXal 

On  substituting  this  into  in.6.6,  we  obtain 

iTCj)  =(1-  t  )’Mpi)  =  Z  (olpi>la)  ni.7.9 


The  coherent  scattering,  in.6.19,  is  then  given  by 


I  Jt'T  )  (  Jl  j  I 


Ijit)  {nj\ 


V  (^1  Pi) 
oo'  1-NT^ 


(q’  I  Pi)* 
1  -  NT*^. 


III.7.10 


We  will  also  evaluate  the  incoherent  scattering  contribution  to  first  order  in  the 
generating  operator  L  .  From  ni.6.18  in  the  limit  of  small  L  ,  we  have  to  first  order 

-  -  ni.7.11 

I  TCt  >  (  tct  I  =  (  1  +  )  I  • 

We  show  in  appendix  A  that  the  first  order  incoherent  term  is  given  by 


I  Tlj  )  (  Jtj  I  - 

T  (7’<7"-y.:d2 

0^0"  1  -  NToNT*o' 


I  HX )  {  Jtj  I  =  I<  I  ^ 


<q"  I  Pi> 

l-NT„. 


(a''+CT’-(T  I  Pi) 
1  -  NT*_.  _ 

CT  -KJ  “C 


I  a  >  ( a'  I. 


ni.7.12 


For  evaluating  the  scattering  cross  section,  we  will  be  interested  in  computing  the  matrix 
elements  (7^1 ...  li^  )  of  the  quantities  111.7. 10  and  ni.7.12.  This  will  be  addressed  in 
section  in.8,  below. 


To  conclude  the  current  section,  we  discuss  the  calculation  of  the  matrix  elements  T(j  and 
To<j>.  First  we  compute  T0.  The  definition  of  Tq  is 

<k„l  '^(0)  lk„>_ 


which  can  be  written  in  configuration  space  as 

T0  =  Jd3r  {  k^lrXr  l1'(7)  lk^>^ 


ni.7.13 


where 


<  rik, 


:>-vV 


r 


(kJr)*. 


in.7.14 


The  quantity  ( r  I  T(  0  )  I  k  ^ )  is  the  effect  of  the  single  scatterer  on  a  simple  plane 
wave,  as  expressed  in  the  configuration  space  representation.  This  is  therefore  the  same 
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quantity  as  was  evaluated  as  the  exterior  solution  in  equations  HI. 3. 15  for  the  single 
uniform  sphere  scatterer,  and  is  given  as 


( ■?  I  ^  ■?  )  I  k  Cno  hn^^^Oc^r)  Pn(cos0) . 

On  substituting  III.7.14  and  111.7.15  into  III.7.13,  we  get 

To  =  J  d3r  e'^*^  ■  Cnohn^^)(kor)P„(cos0) 


ni.7.15 


=  ^  Cn„  /  d3rPn(cos0)  c-iVcos0 
where  an  arbitrary  choice  of  the  configuration  space  coordinate  system  orientation  has 

I— 11^  ■  ^ 

been  made,  placing  the  z  axis  parallel  to  ,  so  that  k^-  r  =  k^r  cos0  ,  which 
simplifies  the  computations. 


Next,  we  evaluate  the  off  diagonal  terms,  Once  again  we  start  from  the  definition 


Toa’  -  ^(®) 


=  Jd3r  <irj7><7i^(7) 


ra.7.17 


This  time  we  allign  the  z  axis  along  k  .  so  that 


(r  I  t(  0  )  i  k^,  )  =  c„^ h„(»(k^,.r) P„(cos9) 


in.7.18 


It  is  convenient  to  write  r  and  k^  in  polar  coordinates,  with  reference  to  figure  3.  In 
particular,  we  have 


(cartesian)  ^y’ 


(polar) 


(k<j,  a,  p) 


53 


Figure  3.  Coordinate  System 


On  explicitly  writing  the  cartesian  components  in  terms  of  the  radial  components  and  the 
polar  and  azimuthal  angles,  we  obtain: 


k<j  =  kgSinacosP 
kcy  =  sina  sinp 
=  k<jCosa 


X  =  rsin9cos(|) 
y  =  r  sin0  sin(]) 
z  =  rcosG 


so  that 


k  •  r  =  k<jr  [sina  sin0  cos((j)-P)  +  cosa  cos0] 


=  k^jrcosy,  where  y  is  the  angle  between  r  and 


Therefore  we  finally  obtain 


Too-  =  Jo  ^  Pn(cos0)  hn^^^Cko-r) . 


ni.7.19 


Section  111.8 

Far  Field  Solution  and  Scattering  Cross  Section 


Far  Field  Solution 

The  approach  used  in  this  section  for  characterizing  the  far  field  solution  makes  use  of 
Kirchoff  s  method  for  expressing  a  scalar  field  inside  a  closed  volume  in  terms  of  the 
value  of  the  field  and  its  normal  derivative  on  the  boundary  surfaces  of  the  volume.  In 
the  current  instance,  we  will  consider  the  volume  of  interest  to  be  the  portion  of  space 
outside  of  the  turbulence  region,  and  extending  to  infinity.  For  the  purposes  of 
visualization,  we  might  consider  the  volume  of  interest  to  be  as  shown  in  figure  4,  taken 
to  the  limit  in  which  the  outer  box  extends  to  infinity.  We  will  be  interested  in  describing 
the  scalar  acoustic  field  and  its  normal  derivative  at  the  surface  Sj  of  the  turbulence 
volume,  and  at  the  surface  at  infinity,  S2. 


For  the  given  acoustic  potential  P(lr  ,t)  =  e*‘“^jts(7) ,  with  Kgir)  satisfying 
(V^  +  k^)  tC5(T^)  =  0,  the  free  space  Green's  function  will  satisfy 
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(V2  +  k^)  G(  r  ,■?' )  =  .5{  r.-?'  ) .  ni.g.  1 


Si 


Figure  4.  Far  Field  Volume 


Green’s  theorem  allows  the  characterization  of  the  field  7ts(r )  in  terms  of  a  surface 
integral 


)■  G(  ■?,■?’  ) &'■  r’  itjCr ' )  ]  da',  III.8.2 


where  n'  is  the  inwardly  directed  normal  to  the  surface  S  . 


The  free  space  Green's  function  is  given  by 


|Up 

G(  )  =  ^ .  ni.8.3 

4jtR 
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where  ^  =  ~r  -Ir' ,  and  R  =  I  it  I- 


On  inserting  this  expression  into  III.  8.2,  we  find 


TtsCr )  =  *  “  s  ^  ^  '  )]da’ ,  ni.8.4 


where  S  =  Sj  +  $2 ,  as  is  shown  in  figure  4,  V  is  the  volume  enclosed  by  Sj  and  S2 ,  and 
the  normal  n'  is  directed  into  V  from  Sj  and  $2  • 

In  the  neighborhood  of  the  surface  S2  at  infinity,  the  scattered  wave  will  behave  as  an 
outgoing  wave,  with  the  implication  that 


and 


e^ 

itsC*")  ~  7- 
47ir 


1  ^tis(r) 
TtsC'r)  ^ 


(ik-l)  , 


SO  that  the  integral  in.8.4  vanishes  on  S2  at  least  as  fast  as  1/r .  We  may  therefore  reduce 
in.  8. 4  to  an  integral  over  only: 


^s( 


ikR  •  ^ 

A'-  [^'  iZs(r' )  +  ik(l+j^)  ^  TtsC'r '  )]da’ . 


m.8.5 


The  contributions  to  this  integral  are  derived  from  the  results  for  Kj  of  the  previous 
section,  coupled  with  the  fact  that  Jts  “  %  *  • 


If  the  observation  point  is  far  from  the  scattering  region,  then 


G(  7,7' ) 


47tr 


ni.8.6 


and 
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ikr 

^s(^)  ~r 


ni.8.7 


In  addition,  we  have 


,ikr 


^s(r)  =  ^  -  K  e-*''  *’  [-it  A'ji3(T^')-fi'  V’  ;t3(r')]da'.  ni.8.8 

^  47t 


from  which  we  can  identify  that 


F(t  .l4c )  =  -  >  c  ^  [  -ii?' «’  >^(  f ' )  -  )]da' 

471  ^1 


—  K  ■  *'  n'-  [  -ii?  Ks(r' )  -  V  7Cs(7'  )]da' 
471 


ni.8.9 


Differential  Scattering  Cross  Section 


The  differential  scattering  cross  section  is  defined  by 


dc(  k  .kjj^c )  _  power  scattered  in  direction  k 


dQ 


ni.8.10 


solid  angle  x  incident  flux  in  direction  kjnc 


where  we  are  interested  in  time  averaged,  rather  than  instantaneous  values  of  flux.  The 
power  scattered  in  direction  is  the  scalar  product  of  the  directed  flux  ^  with  the 
oriented  area  ^  : 


P  =  ^ 


ni.8.11 


The  time  averaged  directed  flux  for  harmonic  fields  is  given  by  (see  Tatarski,  ch  5) 

=  ni.8.12 

where  p  is  the  density  of  air.  We  show  the  selected  coordinate  system  geometry  for  the 
scattering  process  in  figure  5. 
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Figures.  Scattering  Geometry 


We  show  a  plane  wave  incident  on  the  system  from  the  bottom,  with  wave  vector  parallel 
to  the  z  axis,  and  an  outgoing  wave  incident  on  a  small  area  element  on  a  spherical  shell, 
a  distance  R  from  the  center  of  the  scattering  volume.  The  required  factors  are: 


Area  element:  da 
Power:  P 

field  gradient: 


R2dnf , 

R2  dfl  ...  (  //  r  ) 

,  pikr  ,  . 

V  F(t,ri,e) 


=  F(t.ri^)(ik  Y-- 


*ikr 


=  ik  F(  k  ,kinc )  -(at  large  r ,  such  as  r  -  R, ) 
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On  substituting  these  values  into  III.8.12,  we  obtain  the  time  averaged  power  per  unit 
area  scattered  in  direction  Tc 


,-ikR  ^  ,ikR 

V) 


-  ^  t  F*F  DI.8.13 

2R2 


Similarly  we  have  for  the  incident  flux 


(STnc)  =  ^  Im(e’*incZ  il^c  e^inc^) 


0^  p  in.8.14 

Kinr 


On  substituting  these  values  into  111.8.10,  we  obtain 


da(  k  .kjnc )  ^  p*p  in.8.15 

dn 

where  F  =  F(Ti  )  is  the  scattering  amplitude. 

Summary  of  Results  for  Scattering  Cross  Section 

The  final  result  of  the  scattering  cross  section  calculation  is  the  relation  in.8.15 

da(t  .ktnc)  ^ 
dn 

where  F  =  F(lc  ,k^c )  is  the  scattering  amplitude.  The  remaining  task  is  to  substitute 
the  prior  results  derived  from  the  self  consistent  Held  approximation,  and  evaluate  these 
results.  In  particular,  we  have  for  the  scattering  amplitude,  the  expression  in.8.9 


60 


F(t  .l4c )  =  —  K  ^  ^'*1  "sCr ' )  -  JtsCr  '  )]da’ , 
4jt 


from  which  we  can  deduce  that 


F*F 


1  f  f  .-il?- Ji?- 

■(4k)2^s/Si® 


X  A’-  [-ii?  -  ]jis(r  ’ )  [i^  -  ]7ts*(  r  ” )  da’  da" 


— /c  Jc  e*'^‘  «'•  ]  ^s(^'  )^s*Cr  " )  da’  da" 

(47t)2  Si  Si 


=  -irJ,L 

(47t)2  Si  Si 


+  iA'.  t  A"- 


.  A|  A|| 

in*-  V  n**-  k 


+  A*.  A”.  ] 


X  %s(r*  )  da’  da" , 

m.8.16 


where  we  have  interchanged  the  order  of  differentiation  and  averaging,  since  these  refer 
to  different  variables. 
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The  averaged  product  of  the  scattered  waves  can  be  written  as 


JtsCr  '  " )  =  <'r '  I  (  I  JCt  >  -  •  Pi )  )  ( <  •  -  <  Pi  I )  I  ^"  > 


=  {‘r’  iTtxXtCxl  ■?'•> 

-  <r'  i^XPil  r"> 

-  (1^'  IPi><^  r") 
+  ("r '  I  Pi  X  Pi  I  > 


ni.8.17 


From  IIL7.1 1,  we  have  to  first  order 


C?'!  lltT)(ltTll<^">  =  <^'l(l+  1  )lltT>(<'Tl  ■?"> 


=  (■?'  iiiXitT'  ■?">  +  <!?' 11  IhtX^tI  ■?"> 
From  in.7.9,  we  have 

<■?'  =  I  .  (alpiX'r '  la) 

and  similarly, 

<i;i  t")  =  (z  Yrm-  (aiPi><r"io>)*. 

a  1  ^  a 

From  ni.7.12,  we  have 


Cr’IL  IhtX^^tI  "r") 


=  N 


<?a'  a" 
GOO 


'roo"T  0  0"^^'^ 

1  -  NT<,Nr‘o. 


(<y"  I  Pi) 
1-NT^» 


{a"+a'-a  I  pj) 
l-NT* 


Cr' I oXa’ !■?••) 
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The  quantity  <  r  '  I  a  >  is  given  by 


and  similarly, 


The  quantity  <  al  p,)  is  given  as 

(alpi)=  iydh  (  k^\7){7\ 

=  y  J^d^r  *■  e**'inc* 

=  V  lyd^r  e**(*^o  **^inc)-  *■ 

An  issue  of  importance  is  whether  the  incoming  wave  can  be  treated  as  a  normal,  mode  of 
the  system.  This  will  depend  somewhat  on  the  partcular  situation.  If  the  wavelength  of 
the  incident  wave  is  small  compared  to  the  total  size  of  the  turbulence  distribution,  then 
the  density  of  modes  in  the  vicinity  of  the  incident  wavevector  will  be  high,  and 
assumption  that  he  incoming  wave  is  a  normal  mode  should  give  reasonable  results.  This 
condition  is  likely  to  be  satisfied  for  the  application  of  interest,  for  which  case  we  may 
use  the  relationship 


=  5( 


). 


which  significantly  simplifies  the  analysis. 


The  matrix  elements  T^  and  T^o’  ffomin.7.16  andin.7.19,  are  given  as 

To=  ^  Sq  Cn^  J  d3r  h„^^>(k<,r)  PnCcosO)  , 

and 

Too-  =  nS)  ^  e'^Vcosy  hn<*>(koT)  Pn(cose) . 
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Section  111.9 


Relationship  to  Tatarski's  Approach 


In  this  section  we  show  how  Tatarski's  result  for  the  average  value  of  the  flux  density 
vector  of  the  scattered  energy  (averaged  over  one  period  of  oscillation)  follows  from  our 
more  general  formulation  by  keeping  only  lowest  order  terms  in  the  scattering  potential. 
As  in  the  first  part  of  section  III,  we  neglect  the  motion  of  the  medium,  for  simplicity. 
Tatarski's  equations  (5.5),  (5.8),  (5.10),  and  (5.13)  are  then 


(5.5) 


(V2  +  k2)7C  = 
S  v'jl. 


A  series  solution,  Jt  =  ttg  +  7Ci  +  ... ,  is  then  sought  where 


(5.8) 


=  AqC 


iic  r 


(5.10) 


1  ..ikr  ,  .1  ^  ^ 


d^r 


(far  field) 


The  average  scattered  energy  flux  is  given  by 


(5.13) 


S'  =  ^  Im  (7ti*Vjii) 


In  our  more  abstract  formalism,  these  same  equations  are  written  as 


(5.5) 

=  O'  lit) 

ni.9.1 

(5.8) 

ItIq  )  ;  with  <  r'  iTtg  )  =  Aq 

(5.10) 

=  toO’  IJto>  (far  field) 

m.9.2 

(5.13) 

in.9.3 
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As  is  evident  from  equation  (5.10),  or  equivalently  relation  in.9.2,  the  solution  for  k  is 

found  to  first  order  in  the  perturbing  potential  v',  and  the  flux  vector  ^  is  therefore 
obtained  in  the  Bom  approximation.  We  will  show  that  the  relation  III.9.3  for  the 
average  scattered  energy  flux  follows  from  the  results  of  our  choice  of  multiple  scattering 
formalism  in  the  limit  of  neglecting  the  coherent  scattering  entirely  and  keeping  the 
incoherent  results  only  to  lowest  order.  In  particular',  our  results  from  section  III.6  are 


and 


ItitXtctI  =(1+I  +L2  +  L3  +...)  I jct >  { TCx I 


=  (  1  -  L  )*^  (tijI, 


ni.9.4 


\Kj)  =  (  1  +  t  +  (  t  )2  +  ...  )  Ipi) 

=(i-  T )-hpi). 


ni.9.5 


where 


L  =  (  1  -  R  )’^  (M  )  (for  uncorrelated  scaterers) 
=  (1  +R  +R  2  +R  3  +  ...)M 


ni.9.6 


and 


R  X  s  t  X  t  ^ 


m.9.7 


M  X 

and 


=  N  J  d3bj  piC^j)  tC^j)X  t  ^C^j)  =  Ntxt  , 


m.9.8 


ts  =  gov'd- loO-)*! 

=  goO-+toO'JoO-+toO’toO-|oO'+, 


m.9.9 
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To  get  to  Tatarski's  first  Bom  approximation  result,  we  take  the  limit  ofv'  going  to  zero, 

keeping  only  the  lowest  power  non  vanishing  terms  in  each  expansion.  In  particular,  this 
means  that  we  take  t  to  be 

t_  A  A, 

=  goV  .  ni.9.10 

In  addition,  we  find  in.9.5  is  replaced  by 

|JCX>  =  IPi). 

so  that 

I’iiT)  wi  2  iPiXPil  . 

We  may  now  evaluate 

+  IJtx)  (^xl» 

where  we  must  now  consider  I  to  be  given  by  I  s  Af  ,  where 

MX  =  . 

If  we  use  the  result  m.  9. 11,  we  find 

ihjXhti  =  iPi>(Pii+  to  O'  ipiXpil  O'ty  . 

If  we  now  consider  the  scattered  component  only,  we  find  to  this  order  of  approximation, 

I  JCs  )  ( •  =  I  Jtx  )<  I  -  IJCx)  (  Pi  I  -  I  Pi  >  <  TCx  •  +  I  Pi  ><  Pi  I 

=  goO*  ipiXpil  O'ty  ^  in.9.13 

which  is  in  agreement  with  Tatarski’s  result 
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Sction  ni.lO 


Numerical  Algorithm 

The  numerical  algorithm  which  is  presented  in  this  section  is  a  fully  operational  code  for 
obtaining  the  differential  scattering  cross  section  in  the  mean  self  consistent  field 
approximation  for  a  collection  of  randomly  positioned  eddies  confined  to  a  cubical 
turbulence  volume.  The  code  follows  the  theoretical  treatment  developed  in  sections 
in.1.3  through  III.1.8  of  this  document.  Both  the  computation  of  the  coherent  field  and 
the  incoherent  field  have  been  coded.  The  incoherent  field  portion  would  lend  itself  quite 
well  to  operation  as  a  vectorized  code.  It  does,  however,  involve  an  extra  three 
dimensional  momentum  space  integration,  and  operates  very  slowly  on  the  pc  computer 
which  has  been  used  for  code  development .  For  this  reason,  the  incoherent  field  portion 
of  the  algorithm  has  been  commented  out  in  the  phase  I  dgorithm.  This  is  clearly  marked 
by  the  use  of  a  "c*"  notation  in  the  first  column  of  the  listing.  It  could  easily  be  reinstated 
for  evaluation  on  a  vector  processor,  if  desired. 

The  treatment  in  the  numerical  algorithm  makes  use  of  a  code  segment  obtained  from  a 
text  by  Barber  and  Hill,  for  computation  of  the  Bessel  function  expansions  identified  in 
the  T-matrix  plane  wave  representation  portion  of  part  III  of  this  report.  It  does,  however, 
diverge  from  the  coordinate  system  choice  identified  in  section  III.7  of  the  theory.  This  is 
because  an  election  was  made  to  treat  a  definite  cubical  turbulence  volume,  for  which 
case  there  is  a  preferred  orientation  of  the  spacial  coordinates  along  the  axes  of  the  cube, 
rather  than  choosing  the  z  axis  aligned  to  the  k-space  wave  vector.  The  angle  between 
the  space  vector  and  the  wave  vector  is  still  computed,  and  the  result  is  used  in  the 
spherical  expansion  of  the  single  particle  solution.  The  only  signigicant  consequence  of 
the  coordinate  system  orientation  choice  is  that  the  volume  integrals  are  earned  out  in 
cartesian,  rather  than  polar,  coordinates. 

The  selection  of  a  cube  for  the  shape  of  the  turbulence  volume  is  arbitrary,  but  does  lend 
itself  to  computational  simplicity.  There  is,  at  this  stage  of  code  development,  no 
specified  a  priori  shape  to  associate  with  the  physical  turbulence  volume.  Shapes  of  a 
rectangular  prism  character,  but  with  non  cubical  aspect  ratio  could  be  computed  with 
negligible  perturbation  to  the  algorithm  as  it  stands;  many  of  the  required  code  lines  for 
this  modification  have  already  been  inserted,  but  variables  in  these  lines  have  been  set  to 
default  cubical  values. 

The  treatment  of  an  arbitrarily  shaped  turbulence  volume  lies  outside  the  scope  of  phase  I 
development,  but  should  be  relatively  straightforward  to  accomplish  under  a  phase  n.  In 
particular,  a  natural  way  to  approach  the  task  would  be  to  approximate  the  volume  with  a 
collection  of  rectangular  prisms,  joined  at  the  boundaries  by  transport  boundary 
conditions.  This  would  allow  preservation  of  the  main  features  of  the  phase  I  code.  The 
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process  of  joining  adjacent  prisms  would  be  very  similar  in  character  to  the  processes 
used  in  finite  element  computational  fluid  dynamics  algorithms. 

Some  compromises  have  been  made  in  the  development  of  the  phase  I  algorithm,  in  the 
interest  of  both  schedule  and  run  time.  In  particular,  the  code  as  presented  treats  only  a 
single  turbule  size.  The  treatment  of  a  collection  of  sizes  is  a  straight  forward 
modification  in  the  T-matrix  calculkations,  in  which  the  averaging  of  the  collection  of 
waves  to  obtain  a  mean  wave  would  yield  not  a  simple  multiplication  factor  of  N,  but 
rather  a  sum  over  sizes  of  the  product  of  N  for  each  size  with  the  quantity  being  averaged. 
The  resulting  T  matrices  would  have  additional  indices  for  the  turbule  size.  An  algorithm 
for  the  computation  of  the  number  densities  for  each  size,  within  the  confines  of  a 
particular  turbulence  model,  has  been  presented  in  Part  II  of  this  report,  but  has  not  yet 
been  integrated  into  the  scattering  algorithm.  We  have  also  developed  a  code  module  for 
inverting  the  matrix  which  results  from  the  part  II  algorithm  when  a  large  number  of  size 
bins  have  been  selected.  This  code  is  also  operational,  but  has  not  yet  been  fully  tested, 
and  is  not  presented  in  this  phase  I  report  It  is  suggested  that  the  phase  II  code 
development  should  include  as  a  first  step  the  integration  of  these  two  codes  with  the 
body  of  the  algorithm,  and  the  development  of  the  T-matrices  when  there  is  more  than 
one  size. 

Use  of  Algorithm  for  Shadow  Zone  Analysis 

Several  options  present  themselves  for  the  intended  use  of  this  code  in  shadow  zone 
analysis.  One  of  these  options  is  to  use  the  algorithm  as  is,  but  embed  it  in  a  larger  code 
which  has  been  designed  to  simulate  the  shadow  zone  environment  Some  characteristics 
of  the  code  should  be  identified  in  the  assessment  of  the  viability  of  this  option.  First,  the 
election  to  use  a  cubical  turbulence  volume  will  lead  to  volume  boundary  artifacts  when 
the  scatterer  density  is  high,  because  the  planar  surfaces  will  behave  as  mirrors.  This  is 
probably  not  a  significant  issue  at  low  eddy  density,  because  in  that  case,  the  acoustic 
wave  would  penetrate  well  into  the  volume,  and  the  scattered  wave  would  be  principally 
influenced  by  volumetric,  rather  than  surface  effects.  Another  consideration  is  that  a  far 
field  approximation  has  been  used.  This  will  be  highly  appropriate  when  the  turbulence 
volume  is  small  and  reasonably  distant  If  the  smallness  condition  is  not  satisfied,  but 
the  turbulence  volume  is  still  reasoonably  distant  it  may  be  approprate  to  subdivide  the 
volume  into  smaller  regions,  and  compute  the  contribution  from  each  within  the  far  field 
approximation.  This  procedure  will  be  invalid  for  wavelengths  which  are  comparable  to 
the  size  of  the  total  turbulence  volume,  but  such  wavelengths  arc  not  likely  to 
characterize  the  acoustic  spectrum  of  interest  ,  since  this  lies  in  the  midrange  audio 
spectrum. 

A  second  option  in  the  use  of  the  code  in  shadow  zone  analysis  would  be  to  conduct  the 
additional  development  required  to  handle  the  near  field.  This  would  allow 
characterization  of  scatter  from  turbule  distributions  with  a  large  angular  subtense. 
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Additional  theoretical  development  would  be  required  before  a  code  module  of  this 
character  could  be  written. 

Some  consideration  should  be  given  as  to  how  to  interface  the  current  module  with  a 
system  level  environmental  description  model.  One  method  for  treating  the  system  level 
model  would  be  essentially  a  finite  element  fluid  model  with  acoustic  transport.  If  this 
style  of  code  is  envisioned,  the  use  of  the  current  module  would  be  straightforward  in  that 
one  or  more  of  the  finite  elements  would  simply  be  replaced  by  the  current  algorithm  in 
so  far  as  that  elements  acoustic  properties  are  concerned.  If  the  system  level  algorithm  is 
to  be  a  more  simplified  description,  such  a  s  a  linear  circuit  model,  some  thought  would 
need  to  be  applied  as  to  how  to  represent  this  current  algorithm  as  an  equivalent  circuit 
element. 

Finally,  if  the  goal  is  to  obtain  as  literal  as  possible  an  environmental  description,  it 
would  be  desirable  to  extend  the  current  algorithm  to  arbitrary  shapes  of  the  total 
turbulence  volume.  This  would  involve  analysis  of  interface  acoustic  transport  equations, 
and  the  encoding  of  these  equations  between  rectangular  elements.  A  Finite  Volume 
Elements  approach  would  lend  itself  well  to  this  task.  As  an  alternative,  it  might  be 
desired  to  approximate  the  actual  physical  turbulence  volume  by  some  relatively  simple 
three  dimensional  shape  of  high  symmetry,  such  as  a  truncated  cylinder,  an  oblate 
spheroid,  etc.  In  this  latter  case,  some  additional  development  woric  would  be  required  to 
recast  the  integrals  into  the  selected  geometry. 

Phase  I  algorithm 

The  text  of  the  phase  I  algorithm  is  presented  below. 


program  acstc 

complex  hankel(200),hprime(200),xi, 

+  b(200),c(200),e(200),f(200),cn0(200),dum 

common  /avkncom/  N,Tix,Tiy,den 
common  /tcom/  el,nmax,pijmax,cn0 
dimension  xk(3),xld(3),a(200),d(200) 
c*  dimension  3isg(3),isg(3) 
c  N  is  the  number  of  eddies  in  the  volume 
c  TiXjTiy  are  the  real  and  imaginary  parts  of  the  T  matrix 
c  den  is  an  intermediate  expression  used  in  evaluating  T 
c  el  is  the  side  length  of  the  (assumed  cubical)  turbulence  volume 
c  nmax  is  the  highest  order  to  use  in  spherical  expansions 
c  jmax  is  the  number  of  cells  per  side,  when  dividing  up  the 
c  volume  for  doing  volume  integrals, 
pi  =  3.141592653 

c  hankel(n)  is  a  complex  array  for  storing  hankel  functions 
c 
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c  open  input  data  file 

open(  1  ,file='acstc .  inp') 

read(l,*)  maxsg.el,(xk(i),i=l,3),(xki(j)o=l,3),N,jmax,nmax 
c 

c  open  output  file 

open(7,fiIe=’acstc.out') 

c 

c  specify  the  range  of  allowed  values  of  sigma  to  compute 
maxsgl  =  maxsg 
maxsg2  =  maxsg 
maxsgS  =  maxsg 
c 

c  evaluate  the  coefficients  to  be  used  in  spherical  expansions 
c 

c  ...first  some  assumptions,  to  use  as  a  definite  example... 
c 

xa  =  10. 

c  xa  is  the  radius  of  turbule  (cm) 
c 

xldmag=sqrt(xki(  1 )  *xki(  1  )+xki(2)  *xki(2)+xki(3)*xki(3)) 
c  xkimag  is  the  scalar  magnitude  of  the  incident  wave  vector 
c 

rho  =  .0018 

c  rho  is  the  density  of  air  (g/cm'^3)  outside  eddy 
c 

rhoprm  =  rho*  1.0001 
c  rhoprm  is  the  density  inside  the  eddy 
c  (assume  a  small  (0.01  percent)  density  fluctuation) 
c 

xkprm  =  1.0002*xkimag 

c  xkprm  is  the  magnitude  of  the  wave  vector  inside  the  turbule 
c  (assume  a  small  (0.02  percent)  refractive  index  fluctuation) 
c 

c  note:  variations  in  rho  and  n  may  or  may  not  be  related, 
c  depending  on  temperature,  spe^,  moisture,  etc  considerations 
c 

c  The  coefficients  cnO  can  be  found  by  solving  the  relations 
c  a*an0  -  b*cn0  =  c 

c  and  d*an0  -  e*cn0  =  f 
c  with  the  coefficients  a  through  f  given  below 

xi  =  (0.,l.) 

c  xi  is  the  imaginary  number  square  root  of  *1. 
c 

c  evaluate  the  required  (exterior)  bessel  functions  at  the  boundary 
write(6,*) '  start  1st  bcsh  call' 
write(7,*) '  start  1st  besh  call' 
call  besh(xkimag*xa,hankel,nmax+l) 

evaluate  the  required  derivatives  of  bessel  functions  at  boundary 
write  (6,*) '  return  from  besh...  start  derivatives' 
write  (7,*) '  return  from  besh...  start  derivatives' 
hprime(l)  =  -  hankel(2) 
do  40  i=2,nmax+l 
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40  hprime(i)  =  hankel(i-l)  -  (i+l)*hankel(i)/xkimag/xa 
c 

c  write  the  exterior  coefficients  in  terms  of  bessels  and  derivatives 
write(6,*) '  done  with  derivatives,  start  exterior  coefficients' 
write(7,*) '  done  with  derivatives,  start  exterior  coefficients' 
do  50  i=l,nmax 
ip  =  i+1 

b(i)  =  rho*hankel(ip) 

c(i)  =  rho*(l  .+2.*i)*xi**i*real(hankel(ip)) 
e(i)  =  xkimag*hprime(ip) 

50  f(i)  =  xkimag*(l.+2.*i)*xi**i*real(hprime(ip)) 
c 

c  evaluate  the  required  (interior)  bessel  functions  at  boundary 
write(6,*)  ’  done  with  exterior...do  interior  besh  call' 
write(7,*) '  done  with  exterior...do  interior  besh  call' 
call  besh(xkprm*xa,hankel,nmax+l) 
c 

c  evaluate  the  required  derivatives  of  bessel  functions  at  boundary 
write(6,*)  ’  done  with  besh,  interior  start  derivatives' 
write(7,*)  ’  done  with  besh,  interior  start  derivatives’ 
hprime(l)  =  -  hankel(2) 
do  60  i=2,nmax+l 

60  hprime(i)  =  hankel(i-l)  -  (i+l)*hankel(i)/xkprm/xa 
c 

c  write  the  interior  coefficients  in  terms  of  bessels  and  derivatives 
write(6,*) '  evaluate  interior  coefficients  and  solve  for  cnO' 
write(7,*)  ’  evaluate  interior  coefficients  and  solve  for  cnO' 
do  70  i=l,nmax 
ip  =  i+1 

a(i)  =  ihoprm*real(hankel(ip)) 
d(i)  =  xkprm*real(hprime(ip)) 
c 

c  solve  for  the  cnO 

dum  =  b(i)*d(i)-e(i)*a(i) 
if(dum.eq.(0.,0.))go  to  99 
70  cn0(i)=(f(i)*a(i)-c(i)*d(i))/dum 

c 

c  compute  first  term  of  dsigma/dorr«ga 

write(6,*) '  done  with  cnO...  start  dsigma/domega  #1' 
write(6,*)  'cnO  =  ',(cn0(i),i=l,nmax) 
write(7,*)  'cnO  =  ',(cn0(i),i=  1  ,nmax) 
write(7,*)  ’  done  with  cnO...  start  dsigma/domega  #1' 
fil  =  el*el*el*avkna(xki)*fi(xk,xki)/16./pi/pi 
write(6,*)  'fil=',fil 
write(7,*)  'fil=’,fil 
c 

cc*  compute  second  term  of  dsigma/domega 
cc*  sum  fi2  over  ksigma  for  allowed  values  of  sigma 
c*  write(6,*) '  now  do  dsigma/domega  #2' 
c*  write(7,*) '  now  do  dsigma/domega  #2' 
fi2  =  0. 

c*  do  300  il  =  -maxsgl,  maxsgl 
c*  isg(l)  =  il 

c*  do  200  i2  =  -maxsg2,  maxsg2 
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c*  isg(2)  =  i2 

c*  do  100  i3  =  -maxsg3,  maxsg3 

c*  isg(3)  =  i3 

cc*  compute  ksigma 
c*  do400i=l,3 

c*400  xksg(i)  =  2*pi*isg(i)/el 

c*  fil  =  ri2+el*el*el*avkns(xksg,xki)*fi(xk,x]csg)/l  6/pi/pi 

c*  write(6,*) '  done  with  avlms  for  ksigma  =',xksg 

c*  wTite(7,*)  ’  done  with  avkns  for  ksigma  =’,xksg 

c*100  continue 

c*200  continue 

c*300  continue 

c 

c  combine  first  and  second  term  of  dsigma/domega 
c*  write(6,*)  ’  add  fil  and  fi2  to  get  total  dsigma/domega’ 

c*  write(7,*)  ’  add  fil  and  fi2  to  get  total  dsigma/domega' 

fstrf  =  fil +fi2 
write(6,*)  'fstrf  =  '.fstrf 
write(7,*)  'fstrf  =  fstrf 
c 

c  all  done, 
stop 

99  write(6,*) '  boundary  matrix  is  singular' 
stop 
end 


function  fi(xkl,xk2) 

dimension  xkl(3),xk2(3),smallq(3),bigQ(3) 

write(6,*) '  enter  function  fi' 

xf=  1. 

qdotQ  =  0. 

do  100i=l,3 

smallq(i)  =  xkl(i)  -  xk2(i) 

bigQ  =  xkl(i)  +  xk2(i) 

qdotQ  =  qdotQ  +  smallq(i)  *  bigQ(i) 

xf  =  xf*xg(smallq(i)) 

100  continue 

fi  =  qdotQ*qdotQ*xf 
write(6,*) '  done  with  function  fi' 
return 
end 


function  xg(x) 
write(6,*) '  start  function  xg' 
if(x.le.0.001)go  to  10 
xg  =  sin(x)/x*sin(x)/x 
write(6,*) '  done  with  xg' 
return 

10  xg  =  (l.-x*x/6.)*(l.-x*x/6.) 
write(6,*) '  done  with  xg' 
rettim 
end 


function  avkna(xki) 


72 


common  /avkncom/  N,Tix,Tiy,den 
dimension  xki(3) 

write(6,*) '  start  function  avkna','  xki=',xki 
call  T(xki,Tix,Tiy) 

den  =  (1.-  N*Tix)*(l.-  N*Tix)  +  N*N*Tiy*Tiy 

avkna  =  N*N*(Tix*Tix  +  Tiy*Tiy)/den 

write(6,*) '  done  with  avkna* 

return 

end 


function  avkns(xksg,xki) 

common  /avkncom/  N,Tix,Tiy,den 

dimension  xksg(3),xki(3) 

write(6,*) '  start  function  avkns’ 

write(6,*) '  xki  =  *,xki 

writeC?,*) '  xki  =  ',xki 

call  T(xksg,Tsgx,Tsgy) 

write(7,*) '  Tsgx  =  ',Tsgx,Tsgy  =  ',Tsgy 

den2  =  1.-  N*N*(Tsgx*Tsgx  +  Tsgy*Tsgy) 

call  T2(xksg,xki,T2sgx,T2sgy) 

avkns  =  N*(T2sgx*'nsgx  +  T2sgy*T2sgy)/den2/den 

write(6,*) '  done  with  avkns* 

return 

end 


subroutine  T(xk,Tx,Ty) 
common  /Tcom/  el,nmax,pi  jmax,cnO 
complex  sphsum,hn(200),Texp,TT,cnO(200) 
dimension  xk(3),pn(200) 
c  xk  is  the  wavevector 
c  Tx  is  the  real  part  of  TT 
c  Ty  is  the  imaginary  part  of  TT 
c  el  is  the  side  length  of  the  cubical  turbulence  volume 
c  nmax  is  the  cut  off  for  the  summation  of  spherical  functions 
c  jmax  is  the  #  of  cells  into  which  each  side  of  the  cube  is  divided 
write(6,*)  ’  enter  subroutine  T,*  xk=*,xk 
write(6,*)  ’  establish  limits  of  r  integration* 
nxmax  =jmax 
nymax  =  jmax 
nzmax=jmax 
write(6,*)  ’  set  d-cubed  r* 
deltax  =  el/nxmax 
deltay  =  el/nymax 
deltaz  =  el/nzmax 
write(6,*)  '  set  starting  r* 
xO  =  (nxmax  +  l)*deltax/2. 
yO  =  (nymax  +  l)*deltay/2. 
zO  =  (nzmax  +  l)*deltaz/2. 
write(6,*)  ’  find  magnitude  of  k  vector* 
xkmag  =  sqrt(xk(l)*xk(l)  +  xk(2)*xk(2)  +  xk(3)*xk(3)) 

TT  =  (0.,0.) 

write(6,*)  'start  r  integration* 
do  400  iz=l, nzmax 
z  =  iz*deltaz  -  zO 
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Z^Z 

write(6,*)  ’  z=',z,'  z2=',z2 
do  300  iy=l,nymax 
y  =  iy*deltay  -  yO 
y2  =  y*y 

write(6,*) '  y=’,y,'  y2=’,y2 
do  200  ix=l,nxmax 
X  =  ix*deltax  -  xO 
x2  =  x*x 

write(6,*) '  x=’,x,'  x2=’,x2 
r  =  sqrt(x2  +  y2  +  z2) 
writc(6,*)  *  n=*,r 

rdotk  =  (x*xk(l)  +  y*xk(2)  +  z*xk(3)) 

write(6,*) '  rdotk=', rdotk 

rk  =  r*xkrnag 

write(6,*) '  rk=',rk 

cosang  =  rdotk/rk 

write(6,*) '  cosang=', cosang 

c  cosang  is  the  cosine  of  the  angle  between  r  and  k  vectors 
sphsum  =  (0.,0.) 
write(6,*) '  doing  r’4x,iy,iz 
write(6,*) '  start  subroutine  p' 
call  p(cosang,pn,nmax) 

write(6,*) '  done  with  subroutine  p,  enter  subroutine  besh' 
call  besh(rk,hn,nmax+l) 

write(6,*)  ’  done  with  besh...  start  sphsum  integral' 
do  100  n=l,nmax 

100  sphsum  =  sphsum  +  cn0(n)*hn(n+l)*pn(n) 

Texp  =  cmplx(cos(rdotk),-sin(rdotk)) 

200  TT  =  TT  +  Texp*sphsum*deltax*deltay*deltaz 
300  continue 
400  continue 

TT  =  Tr/el**1.5 
Tx  =  real(TT) 

Ty  =  aimagCIT) 

write(6,*) '  done  with  subroutine  T' 

return 

end 

subroutine  T2(xkl,xk2,Tx,Ty) 
common  /Tcom/  el,nmax,pi,jmax,cn0 
dimension  xkl(3),xk2(3),pn(200) 
complex  sphsum, hn(200), Texp, T,cn0(200) 
c  xkl  is  ^e  "outgoing"  wavevector 
c  xk2  is  the  "incoming"  wavevector 
c  Tx  is  the  real  part  of  T2 
c  Ty  is  the  imaginary  part  of  T2 
c  el  is  the  side  length  of  the  cubical  turbulence  volume 
c  nmax  is  the  cut  off  for  the  summation  of  spherical  functions 
c  jmax  is  the  #  of  cells  into  which  each  side  of  the  cube  is  divided 
nxmax  =  jmax 
nymax  =  jmax 
nzmax  =  jmax 
deltax  =  el/nxmax 
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deltay  =  el/nymax 
deltaz  =  el/nzmax 
xO  =  (nxmax  +  l)*deltax/2. 
yO  =  (nymax  +  l)*deltay/2. 
zO  =  (nzmax  +  l)*deltaz/2. 

xklmag=sqrt(xkl(l)*xkl(l)+xkl(2)*xkl(2)+xkl(3)*xkl(3)) 

xk2mag=sqrt(xk2(l)*xk2(l)+xk2(2)*xk2(2)+xk2(3)*xk2(3)) 

T=(0..0.) 

do  400  iz=l, nzmax 
z  =  iz*deltaz  -  zO 
z2  =  z*z 

do  300  iy=l, nymax 
y  =  iy*deltay  -  yO 
y2  =  y*y 

do  200  ix=l, nxmax 
X  =  ix*deltax  -  xO 
x2  =  x*x 

r  =  sqrt(x2  +  y2  +  z2) 

rdotkl  =  (x*xkl(l)  +  y*xkl(2)  +  z*xkl(3)) 

rdotk2  =  (x*xk2(l)  +  y*xk2(2)  +  z*xk2(3)) 

rkl  =  r*xklmag 

rk2  =  r*xk2mag 

cosgam  =  rdotkl/rkl 

costht  =  idotk2/rk2 

c  costht  is  the  cosine  of  the  angle  between  r  and  k2  vectors 
c  cosgam  is  the  cosine  of  the  angle  between  r  and  kl  vectors 
sphsum  =  (0.,0.) 
call  p(costht,pn,nmax) 
call  besh(rk2,hn,nmax+l) 
do  100  n=l,nmax 

100  sphsum  =  sphsum  +  cn0(n)*hn(n+l)*pn(n) 

Texp  =  cmplx(cos(rdotkl),-sin(rdotkl)) 

200  T  =Texp*sphsum*deltax*deltay*deltaz 
300  continue 
400  continue 
T  =  T/el**1.5 
Tx  =  real(T) 

Ty  =  aimagCT) 

return 

end 

Subroutine  besh(x,hankel,nc) 
c  This  subroutine  taken  from  Barber  and  Hill,  'Light 
c  Scattering  by  Particles:  Computational  Methods,' 
c  World  Scientific,  vol2, 1990 
c 

c  calculates  hankel  function 
c  bj  =  bessel  function  of  first  kind 
c  by  =  bessel  function  of  second  kind 
c  X  =  real  argument 

c  nc  =  order  of  functions...  starts  at  zero...  ie: 
c  hankel(l)  =  (j0,y0);hankel(2)  =  (jl,yl);etc. 

c 

complex  hankel(nc) 
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dimension  bj(201),by(201),t(3) 
c 

c  by(*)  calculation...  obtain  0th  &  1st  order  functions 
a  =  sin(x)/x 
by(l)  =  -cos(x)/x 
by(2)  =  by(l)/x  -  a 
c 

c  obtain  the  higher  order  functions  by  upward  recursion 
do  n  =  3,nc 
m  =  red(n-2) 

by(n)  =  (2.’*'m+l.)*by(n-l)/x-by(n-2) 
end  do 
c 

c  bj(*)  calculation...  set  the  starting  order  for  downward  recursion 
nst  =  int(x+4.05*x**.3333+2.+(101.+x)**.5) 
c  the  t(*)  array  is  used  to  recur  down  to  the  two  highest  order  functions 
c  that  are  needed, 
c 

c  set  starting  values  for  the  two  highest  orders  nst  and  nst-1 
t(3)=0. 
t(2)=l.e-35 
c 

c  recur  downward  to  obtain  orders  nc-1  and  nc-2 
do  i  =  nst- l,nc- 1,-1 
ri  =  real(i) 

t(l)  =  (2.*ri+l.)*t(2)/x-t(3) 
t(3)  =  t(2) 
t(2)  =  t(l) 
end  do 
c 

c  continue  downward  recursion  to  order  one 
bj(nc)  =  t(3) 
bj(nc-l)=t(2) 
do  i  =  nc-2, 1,-1 
ri  =  real(i) 

bj(i)  =  (2.*ri+l.)*bj(i-i-l)/x-bj(i+2) 
end  do 
c 

c  calculate  the  scale  factor  and  the  functions 
alpha  =  a/bj(l) 
do  k  =  l,nc 

hankel^)  =  cmplx(bj(k)*alpha,by(k)) 
end  do 
return 
end 

subroutine  p(x,pn,nmax) 
dimension  pn(200) 
pn(l)  =  x 

pn(2)  =  0.5*(3.*x*x-l.) 
do  10  i=3,nmax 

10  pn(i)  =  ((2.*i.l.)*x*pn(i-l)  -  (i-l)*pn(i-2))/i 
return 
end 
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Appendix  A 

First  Order  Incoherent  Effects 

We  have  shown  in  section  III.6  that  the  incoherent  radiation  component  is  given  by 

Itij)  1  )*^  • 

with 

L  =(1  -R  )-l(Af  +Q  R  ), 

where  Af  represents  the  contribution  due  to  individual  scatterers  and  Q  -  R  represents 
the  contribution  due  to  two  particle  correlation  effects  among  pairs  of  scatterers.  Within 
the  context  of  our  phase  I  assumptions,  correlation  effects  have  been  deferred  to  the 
future,  so  that  we  have 


L  =(1  -R  , 


A.l 


The  superoperators  R  andM  are  defined  for  any  configuration  independent  operator  X  as 

_  _+  A.l 

R  X  =  t  X  t  ^ 


and 

M  X  =  2jd3rjp(r7)  t(r7)Xt^r7). 


The  averaged  operator  ^  is  given  in  relation  III.6.5  as 

T  =  Zjd3rjp(r7)  t(r7), 

and  has  the  plane  wave  expansion  IIL7.6 

'Y  =  NlTa 


a 


A.4 


The  expression  ( 1  -  R  )*^  M  can  be  expressed  in  the  plane  wave  basis  by  expanding 
the  operator  in  a  series,  and  evaluating  the  separate  terms  of  the  series.  The  result  is  the 
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combined  with  the  plane  wave  expansion  of  M .  The  first  term  in  the  series  is  X  itself, 
which  can  be  represented  in  the  plane  wave  basis  as 

1X1=  EZ  la)  <ai  X  la' )  (o' I  s  Z  Z  X^o- la)  (o' i . 
a  a'  o  a 

The  second  term  in  the  series  isR  X ,  which  can  be  expressed  in  the  plane  save  basis  as 

R  X  s  T  xT*=  Z  Z(NTo)  ioXcSXta'XffI  (NTo*) 

a  o’ 

=  Z  Z  (NTa)  (NTo  *)  X^o'  la)  (a*  I  .  A.5 

o  o' 


The  third  term  in  the  series  is 

K  2  X  =  R  ( R  X ).  A.6 

=  Z  Z  (NTjj)  la)  (ai  ( R  X )  la' )  (o'  i  (NT^  *)  (...from  A.5) 
o  o’ 

=  Z  Z  (NT^j)  la)  (ol  Z  Z  (NTo")  (NT^j-*)  X^-a'-  la")  <a'"  I  a' )  (o'  I  (NT^  *) 
o  o'  o"  o'" 

=  Z  Z  Z  Z  (NTo")  (NTjj.-*)  (NTo)  (NT^  *)  Xo”o-  I  a)  (ala")  (a"'  I  o’ )  (o'  I  A.7 
o  o'  o"  o'" 

=  Z  Z  Z  Z  (NTo-)  (NT^- *)  (NTo)  (NTo  *)  Xo»o-  I  a)  600"  ( a’  I 

o  o'  o"  o’" 

=  Z  Z  (NTo)2  (NTo  *)2  Xoa-  I  a)  ( a’  I .  A.8 


By  induction,  then, 


R  m  X  =  R  ( R  X ) 


=  Z  Z(NTo)  la)(oi(R  ">-1 X )  la’ ) (o' 1  (NT^ •*) 
o  o’ 

=  Z  Z(NTo)"'(NTo->  X^a'  la)(a’l. 
o  o’ 

On  combining  the  terms  in  the  sum  we  get 


A.9 
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( 1  -R  )*^  X  =  ZZ  [l  +  (NTo)(NTo.*)  +  (NTo)2(NTo*)2  +  ... 

+  (NTor(NT<j->  +...  ]x^c' 


a  o 


=  ZZ  - 3—1;-  Xoa*la)(a’l. 

o  o'  1-  (NToXN^o’  ) 


A.  10 


A.ll 


In  order  to  get  I  X  =  ( 1  -  R  )’^  M  X ,  we  need  next  to  compute  M  X  in  the  plane 
wave  basis: 

M  X=  2jd3rjp(r7)t(r7)xf(r7) 

=  NiJd3rj  lt(r7)lXlf(r"*)l. 


The  integrand  of  this  expression  is 
1  t(r7)lXlf(r7)l  = 

=  ZZZZ  loiXail  tCr]")  la2)(o2lX  lajXaal  ) '  <^4X  ^^4 1- 

Gj  O2  G3 

=  Z  Z  Z  Zl  ai>e’^  ’  '‘‘^2  *’j  7^,02  '^*0403  e’*  '‘j(  <^4 ' 

=  Z  Z  Z  Z  e*'  '  ‘^°2  +  ^03  -  *44 )'  *’j  Toi02T*04a3  X02O3  '  <^1>  <  <^4 ' 

Gi  G2  G3  G4 

so  that  the  integral  becomes 
M  X  =  Zj d3rj  p(  r"' )  t (  )  X  )  = 


A.12 


A.13 


A.14 


J 


N 

V 


J  d3nZ  Z  Z  Z  e'"  *’j  Ta-oj  1*0403  X02O3 '  1 


Oj  O2  O3  04 


-  N  Z  Z  Z  Z  5oj+03,O2+O4  '^0102*^0403  X02O3  *  ^lXct4  I 

Oi  O2  O3  O4 

=  N  Z  Z  Z  Toj027  04(02+04-01)  Xo2(02+04-0i)  *  ^lX<^4  I 

Gj  O2  G4 


A.15 


A.16 
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To  continue  the  derivation  of  L  X  =  ( 1  -  R  )'*  M  X  ,  we  substitute  X  ^  M  X  in 
relation  A.IO.  To  accomplish  this,  we  need  to  evaluate  the  matrix  elements  of  M  X 
which  we  derive  from  relation  A.  16: 

(alM  XI(J')  =  NZZ  ZTqjOjT  O4(a2+04-0i)^a2(02+<^4‘<^i)  ^  ^  ^  ^  ) 

Oi  02  04 

=  N  L  Tjj02  T  o'(02+0'-0)  ^02(02+0'-0) 

02 


We  next  substitute  this  value  into  A.IO,  obtaining 

1X  =(»  -R)-'M  X  =SS  ‘  .  (MX)„„.  laXo'l 

o  o  1-  (NT<j)(NT<j'  ) 


=  XIS  - - T 

oo’02  1-(NTo)(NT/) 


NT0CJ2  T  o’(02+0’-0)  ^02(02+0’-0)  ^  ^  * 


A.18 


We  are  intersted  in  evaluating  A.18  when  the  quantity  X  is  chosen  to  be  I  Jt^  >  ( ttj  I . 
The  plane  wave  representation  of  this  quantity  is  given  in  in.7.10.  In  particular,  we  have 


1  Jt-p  )  ( Jtx  ^ 


y  (gjl  Pi) 

ato-  1-NT^^ 


1  -  NT^.. 


laiXa"  I. 


The  required  matrix  elements  are  found  to  be 
( 02  ^1  aj+tf-o)  =  I .  <  02  I  o,  X  o"  1 02+<J'-<J> 

*  Cl  a" 

_  Pi)  (a2+g'-gl  Pi)* 


On  substituting  these  values  into  A.18,  we  obtain 


222 _ i - 

oo'o2l-(NTo)(NTo*) 


MT  T*  ^^2*  Pi) 

tN  I002  *  0'(02+<J’-<J)  1  -  NT 

O2 


{a2-Ky'-al  pj)* 

1  *  NT  (a2+o'-o) 


I  aXo*  I. 


A.20 


Which  is  the  result  we  were  looking  for. 
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